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Abstract 


Classical  Floquet  theory  is  examined  in  order  to 
generate  a  canonical  transformation  to  modal  variables  for 
periodic  systems.  This  transformation  is  considered 
canonical  if  the  periodic  matrix  of  eigenvectors  is 
symplectic  at  the  initial  time.  Approaches  for  symplectic 
normalization  of  the  eigenvectors  had  to  be  examined  for 
each  of  the  different  Poincard  eigenvalue  cases.  Particular 
attention  was  required  in  the  degenerate  case,  which 
depended  on  the  solution  of  a  generalized  eigenvector. 
Transformation  techniques  to  ensure  real  modal  variables  and 
real  periodic  eigenvectors  were  also  needed. 

Periodic  trajectories  in  the  restricted  three-body  case 
were  then  evaluated  using  the  canonical  Floquet  solution. 

The  system  used  for  analyses  is  the  Sun-Jupiter  system. 

This  system  was  especially  useful  since  it  contained  two  of 
the  more  difficult  Poincare  eigenvalue  cases,  the  degenerate 
case  and  the  imaginary  eigenvalue  case.  The  perturbation 
solution  to  the  canonical  modal  variables  was  examined  using 
both  an  expansion  of  the  Hamiltonian  and  using  a 
representation  that  was  considered  exact.  Both  methods 
compared  quite  well  for  small  perturbations  to  the  initial 
condition.  As  expected,  the  expansion  solution  failed  first 
due  to  truncation  after  the  third  order  term  of  the 
expansion . 

xi 


CANONICAL  FLOQUET  PERTURBATION  THEORY 


I.  Introduction 

Analysis  techniques  of  periodic  systems  have  been 
available  for  many  years  since  Floquet's  work  on  time 
periodic  linear  systems  in  1883.  The  primary  use  of  Floquet 
theory  is  to  find  the  characteristic  or  Poincare  exponents 
of  the  system  and  thereby  determine  the  stability  of  the 
periodic  system.  In  more  recent  studies,  Floquet  theory  has 
been  used  to  construct  a  set  of  periodic  modal  vectors  for 
analysis  (Wiesel,  1981;  Calico  and  Wiesel,  1984;  Ross, 

1991) .  Unfortunately,  none  of  these  works  noted  that 
standard  Floquet  theory  is  not  canonical.  In  order  to  find 
a  useful  set  of  periodic  modal  variables,  a  canonical 
version  of  Floquet  theory  must  be  found. 

This  study  will  look  at  the  required  methods  needed  to 
produce  a  canonical  transformation  from  a  time  periodic 
linear  Hamiltonian  system  to  modal  variables  using  Floquet 
theory.  Special  attention  will  be  given  to  the  various 
types  of  Poincard  exponents  encountered  in  these  types  of 
systems . 

One  of  the  many  interesting  applications  for  canonical 
Floquet  theory  is  that  of  periodic  orbits  in  celestial 
mechanics.  As  is  generally  known,  the  two-body  system  "is 
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the  only  gravitational  problem  for  which  a  closed-form 
solution  has  been  found"  (Wiesel,  1989:45).  But  in 
searching  for  exact  two-body  systems  in  our  solar  system 
alone,  it  is  apparent  that  perturbations  due  to  other  bodies 
must  be  considered. 

Therefore,  while  canonical  Floquet  theory  is  then  the 
primary  focus  of  this  study,  a  secondary  emphasis  will  be  on 
the  analysis  of  the  restricted  three-body  orbit.  The 
particular  system  to  be  looked  at  is  the  Sun- Jupiter  system. 
While  by  no  means  the  most  interesting  system  dynamically, 
its  mass  ratio  of  9.5388E-4  makes  the  study  easier  to  handle 
at  this  stage. 
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II.  Historical  Development 


The  initial  groundwork  for  analysis  of  periodic  systems 
was  laid  out  by  Floquet  as  described  in  Chapter  I.  A 
detailed  search  of  the  library's  resources  turned  up  very 
little  information  that  expanded  on  the  work  of  Floquet  or 
would  help  in  this  study.  The  majority  of  the  work  found  in 
the  areas  of  this  study  is  that  of  Dr.  Wiesel  and  Capt  Ross. 

2.1  Canonical  Transformations 

In  Ross  (1991)/  a  method  for  finding  a  set  of  modal 
vectors  from  a  periodic  system  was  defined,  but  the 
canonical  behavior  of  the  transformation  was  not  adequately 
considered.  As  defined,  there  are  actually  two 
transformations  required  to  change  the  original  Hamiltonian 
into  modal  coordinates.  The  second,  or  modal  transformation, 
was  shown  with  great  detail  to  be  canonical  (Wiesel, 
1981:232-234;  Pars,  1965:453-483).  But,  the  first,  or 
Floquet  transformation,  was  not  completely  examined  until 
Dr.  Wiesel' s  discovery  of  a  proof  for  this  transformation 
(Siegel  and  Moser,  1971:99-101).  This  proof  is  presented  in 
Chapter  III.  Armed  with  these  tools,  the  canonical 
transformations  using  Floquet  theory  could  now  be  completely 
tackled. 

2.2  Real  Canonical  Transformations 

In  working  on  this  effort,  it  was  found  that  the 
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standard  transformations  did  not  always  yield  real  modal 
vectors,  and  an  additional  study  in  transforming  complex 
eigensystems  into  real  eigensystems  was  accomplished.  This 
was  a  significantly  easier  task,  in  comparison  to  the  first, 
as  much  has  been  written  in  linear  algebra  and  matrix 
algebra  textbooks  on  these  types  of  transformations. 

2.3  Perturbation  Theory  on  Periodic  Orbits 

The  only  source  for  this  part  of  the  study  was  the  work 
by  Ross  where  he  states  "It  is  unique  to  use  the 
eigenvectors  and  Poincar4  exponents  of  the  periodic 
trajectory,  to  canonically  transform  the  generic  equations 
of  motion  into  nearly-periodic  ones"  (Ross,  1991:3).  While 
Ross's  work  followed  the  theory  of  Dr.  Wiesel,  his  was 
indeed  the  first  to  significantly  analyze  the  possibilities 
of  these  methods. 
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III.  Theory 


3.1  Classical  Floquet  Theory 

In  classical  Floquet  theory,  given  a  Hamiltonian  system 
where  H(X, t)  is  the  Hamiltonian  function,  and  the  state 
vector  is 

. .  .  . ,  q^tpj)  ,  i=l,n  (1) 

The  Hamiltonian  equations  of  motion  are  then  defined  as 

4,.^  and  i-l,n  (2) 


for  each  coordinate,  q^,  and  its  conjugate  momenta,  p^. 
This  can  be  more  compactly  written  as 


where  Z  is  a  correlation  matrix  defined  as 


Z  = 


0  10  0- 

-10  0  0 

0  0  0  1  - 

0  0-10 


(4) 


and  Z  is  of  order  2N  (the  number  of  coordinates  and  momenta 
in  the  state  vector)  (Wiesel,  1981:232,233).  Since  the 
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matrix  Z  is  skew-symmetric  and  has  a  determinant  equal  to 
one,  it  follows  the  well  established  facts  that 
Z*  =  z"!  =  -z.  The  equations  of  variation  then  come  from  the 
linearization  of  the  Hamiltonian  equations  of  motion.  This 
is  accomplished  by  differentiating  Eq(3)  by  the  state  vector 
yielding 


ax  ^ 

■3x  ax' 


A{t) 


(5) 


Here  A(t)  is  the  linearization  of  the  dynamics  and  the 
equations  of  variation  become 

S  =  A(t)X  =  (6) 

ax' 

(ROSS,  1991:10,11),  Here  x  indicates  coordinates  on  the 
tangent  space  of  the  system  defined  by  Eq(3)  (i.e.,  on  a 

trajectory  near  the  periodic  motion) .  From  this,  a  general 
solution,  close  to  the  periodic  system  of  Eq(3),  can  be 
written  using  the  fundamental  matrix  4>(t,to),  as 

X(t)  =  <I>(t,  t„)X(t„)  (7) 

where  the  fundamental  matrix  also  obeys 

^(<D)  =  A{t)<t>  =  (8) 

^  ax' 
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Given  that  the  original  system  does  include  a  periodic 
motion,  then  both  Eqs(6)  and  (8)  are  time  periodic  linear 
differential  equations.  Note  that  the  equations  of 
variation  actually  arise  from  a  first  order  Taylor  series 
expansion  of  the  original  Hamiltonian  about  its  nominal 
trajectory  X(t) 

(9) 

The  main  conclusion  of  classical  Floquet  theory  then  is 
that  the  periodic  fundamental  matrix  4>  can  be  decomposed 
into 

(10) 

where  F(t)  is  periodic  and  J  is  a  Jordan  normal  form. 

The  Floquet  solution  then  is  obtained  by  integration  of 
Eqs(3)  and  (8)  over  one  period  of  the  motion.  One  result  is 
the  monodromy  matrix  <I>(T+to,  t^) ,  where  t  is  the  period.  The 
matrix  F(t)  is  then  the  system  eigenvectors  and  exp{jT}  is 
the  system  eigenvalues  at  one  period.  More  specifically,  J 
is  a  Jordan  normal  matrix  of  Poincard  exponents,  a^.  These 
cOi  then  describe  the  stability  of  the  system  (Pars, 
1965:461-467).  Once  0(T+to,to)  is  found,  and  noting  F(to+t) 

=  F(to)  Eq(lO)  can  be  rewritten 
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4>  =  Fe*^F"^ 


(11) 


where  the  time  indices  on  O  and  F  have  been  dropped  due  to 
their  periodic  nature.  A  standard  math  package  can  now  be 
used  on  4>  to  find  the  system  eigenvalues,  and  the  system 
eigenvectors,  These  eigenvectors  are  placed  in  the 

columns  of  the  F  matrix  and  the  eigenvalues  are  the  diagonal 
elements  of  the  exp{jT}  matrix.  The  diagonal  elements  of  J 
are  then  found  using 

00,  =  lln(X,)  (12) 


or  more  conveniently  for  complex  eigenvalues. 


Re  (00,) 
Jm(co,) 


(13) 


It  is  worth  noting  at  this  time  that  the  Poincar^ 
exponents  occur  only  in  positive/negative  pairs  for  a 
canonical  system.  This  leads  to  exactly  four  possible  types 
of  Hamiltonian  coordinate/momenta  pairs;  1) 
positive/negative  real  oo^,  2)  positive/negative  imaginary 
o\,  3)  a  pair  of  zero  cot,  snd  4)  a  'box'  or 
positive/negative  pair  of  complex  conjugate  cct-  Only  the 
first  three  will  be  examined  in  detail  in  this  study. 
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Completion  of  the  Floquet  solution  at  this  point 
requires  knowledge  of  F(t)  over  one  period.  Differentiating 
Eq(ll)  and  substituting  in  Eq{8)  and  (10)  yields 

Fit)  =  z4!£r(t) -F(t)  J  (14) 

dx" 

A  complete  solution  to  the  equations  of  variation,  through 
Eq(7),  can  now  be  characterized  over  one  period  of  the 
motion.  Knowledge  of  F  over  one  period  also  allows  for  a 
transformation  to  a  set  of  periodic  modal  variables,  §’(t) 
through  the  relationship 

y(t)  =  (15) 

by  the  following  derivation 

Hit)  = 
x(t)  = 

F-Mt)X(t)  = 

y(t)  =  (16) 

To  be  useful  for  further  study,  this  transformation  needs  to 
be  canonical. 

3.2  Canonical  Floquet  Theory 

The  most  important  criterion  for  a  canonical 
transformation  is  that  the  new  Hamiltonian  must  also  follow 
Eq(3),  that  is 
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(17) 


where  K  is  the  new  system  Hamiltonian  and  Y  the  new  system 
state  vector.  In  order  to  ensure  this,  the  original  system 
eigenvector  matrix,  F,  must  satisfy  the  relation 

Z  =  F'^ZF  (18) 


In  other  words,  this  means  that  the  new  Hamiltonian 
equations  of  motion  also  follow  Eq(2)  for  the  new  state 
vector  and  the  same  correlation  matrix.  The  formal  proof 
defines  the  eigenvector  matrix  that  satisfies  Eq(18)  as  the 
symplectic  eigenvector  matrix  (Siegel  and  Moser:  1971,99- 
101) .  As  in  other  work  with  canonical  systems,  this 
transpose  of  F  is  a  standard  transpose  and  not  a  Hermitian 
transpose.  Equation  (18)  for  a  canonical  transformation  is, 
unfortunately,  only  applied  to  constant  eigenvector  systems 
in  Siegel  and  Moser.  Proof  of  its  usefulness  in  a  periodic 
system  must  be  shown  before  this  study  can  continue.  Since 
Z  is  a  constant  matrix, 

i  =  t^ZF+F^Zt  (18) 

must  equal  zero.  Combining  Eqs(5),  (14),  (19),  using  the 

identity  relations  for  Z,  and  assuming  Eq(18)  holds  true, 
(19)  reduces  to 
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-J^Z-ZJ  =  0 


(20) 


Since  J  and  Z  are  constant  over  the  period  of  the  motion, 
direct  substitution  will  prove  Eq(20)  true  for  either  the 
degenerate  or  non-degenerate  case.  Therefore,  Eg (18)  also 
holds  true  for  a  periodic  system  of  eigenvectors  (Wiesel  and 
Pohlen,  1992:7). 

If  F  forms  a  canonical  transformation,  what  is  the  new 
Hamiltonian?  Differentiating  Eq(16)  yields 

y(t)  =  Je-^'^-^^’yCt,)  =  Jyit)  (21) 

where,  as  in  equation  (6) ,  J  must  be  of  the  form  A{t)  or 

J=A(t)  =  z4!5  =  ZS  (22) 

dT' 

With  S  defined  as  d^K/df^,  the  new  equations  of  variation 
come  from  the  new  variational  Hamiltonian 

^  (23) 

(Wiesel  and  Pohlen,  1992:7) 

3.2.1  Symplectic  Normalization  for  Real  Poincare 
Exponent 3 /Independent  Eigenvectors 

For  a  first  look  at  symplectic  normalization,  consider 

the  case  where  each  Poincare  exponent  generates  a  unique 

real  eigenvector.  The  symplectic  eigenvectors,  p^,  can 
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differ  from  the  original  eigenvectors,  by  only  a 
constant  d^.  In  matrix  notation  this  would  be 

E  =  FD  (24) 

with  D  a  diagonal  matrix  of  constant  multipliers,  and  E  the 
symplectic  eigenvector  matrix.  Substituting  FD  into  Eq(18) 
for  F  yields,  after  rearrangement, 

D  iZD-i  =  F^ZF  (25) 


The  form  of  the  left  hand  side  of  Eg (25)  is 


D-^ZD^ 


■3^  ° 


0  0  ••• 

0  0 


0 

0 

0 


1 

0 


i=l,  2N 


(26) 


If  the  di  are  selected  to  satisfy  Eg (25)  then  D  will 
transform  F  into  the  symplectic  S.  This  then  allows  some 
freedom  in  the  selection  of  the  d^.  Siegel  and  Moser 
suggest  the  selection  of  one  for  all  the  odd  d^,  which  then 
uniquely  determines  the  corresponding  even  dj.  This  method 
generates  one  possible  multiplier  matrix,  D,  and  through 
Eq(24)  a  symplectic  transformation  matrix,  E,  for  any  system 
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with  unique  real  (Bj,  (Siegel  and  Moser,  1971:101). 

3.2.2  Symplectic  Normalization  for  Imaginary  Poincar4 
Exponents /Complex  Conjugate  Eigenvectors 

Unfortunately,  not  all  of  a  system's  pairs  can  be 
normalized  as  easily  as  pairs  of  positive/negative  reals 
with  real  eigenvectors.  In  the  case  of  imaginary  pairs, 
where  the  associated  eigenvectors  will  be  complex  conjugate 
pairs,  the  procedure  as  outlined  by  Siegel  and  Moser  would 
result  in  two  different  multipliers  that  would  destroy  the 
complex  conjugate  nature  of  the  eigenvectors.  Fortunately, 
this  problem  is  easy  to  solve.  If  we  look  at  one  particular 
multiplier  pair  entry  in  Eq(26) 

(D-^ZD-^),,  =  ^  =  (F^ZF),,  (27) 


where  the  subscripts  on  r*zr  and  D'^ZD'^  signify  the  row  and 
column  respectively  of  these  matrix  products.  For  this  to 
result  in  a  symplectic  transformation  and  in  order  to 
maintain  the  complex  conjugate  eigenvectors,  d^  must  equal 
d2.  Therefore  the  multiplier  for  both  eigenvectors,  dj^^,  is 
found  simply  from 


d 


1,2 


(  1 
^(F^ZF) 


(28) 


The  fact  that  the  scale  factors  are  coupled  is  not 
surprising.  While  the  eigenvectors  of  a  linear  system 
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can  usually  be  normalized  independently,  in  a  canonical 
problem  they  must  be  normalized  as  canonical  pairs. 
(Wiesel  and  Pohlen,  1992:10) 

Therefore,  this  method  is  also  appropriate  for  use  on  a  pair 
of  eigenvectors  associated  with  real  co^. 

3.2.3  Symplectic  Normalization  for  a  Pair  of  Zero  Poincar4 
Exponents/Repeated  Eigenvectors 

This  study  would  not  be  complete  without  a  look  at  the 

degenerate  case  of  a  pair  of  zero  The  degenerate  case 

commonly  occurs  in  a  Hamiltonian  system  with  conserved 

quantities.  Equally  common  is  a  rank  deficiency  in  the 

eigenvector  matrix  when  a  pair  of  zeros  is  encountered. 

That  is,  there  will  be  a  repeated  eigenvector  (Wiesel  and 

Pohlen,  1992:10). 

3. 2. 3.1  Determination  of  Generalized  Eigenvector 

The  first  task  in  the  symplectic  normalization  of  the 
degenerate  case  is  to  determine  the  generalized  or  extended 
eigenvector  that  removes  the  rank  deficiency  in  F.  The 
generalized  eigenvector,  is  found  using 

=  r,  (29) 


for  a  constant  coefficient  case  (Reid,  1983:346).  But  in  a 
periodic  system  where  the  form  of  exp(jT)  and  the  eigenvalue 
matrix  are,  respectively, 


J  = 


p  1 
0  0 


and 


(30) 
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and  the  more  appropriate  form  of  Eq(29)  is  then  written 


(31) 


As  a  check  on  the  new  F  with  the  extended  eigenvector,  the 
system  should  satisfy  4*'-FA=0.  Note  from  Eq(31)  that  the 
generalized  eigenvector  is  not  entirely  arbitrary  and 
therefore  cannot  be  normalized  in  the  standard  fashion.  It 
will,  however,  be  indeterminable  to  an  additive  multiple  of 
the  repeated  eigenvector.  In  fact, 

Kx  =  (32) 


defines  a  generalized  eigenvector  for  any  value  of  a  (Wiesel 
and  Pohlen,  1992:10). 

3, 2. 3. 2  Normalization  for  the  Degenerate  Case 

Because  of  the  arbitrary  value  of  o,  the  multiplier 
matrix,  D,  is  not  necessarily  diagonal.  Instead  it  will  be 
of  the  form 


D  = 


” 

a 

0  d. 


(33) 


for  the  repeated/generalized  eigenvector  pair.  Since  D  is 
no  longer  diagonal,  Eq(25)  will  also  change  to 

{D^)  =  F^ZF  (34) 
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since  D*  is  no  longer  equal  to  D.  After  the  required  math 
on  the  left  hand  side  of  Eq(34),  the  result  for  a  degenerate 
pair  of  eigenvectors  is  the  same  as  that  for  any  other  pair 
in  Eq  (26)  ,  or 


(D’’) 


1  1 
I 


(35) 


Equation  (35)  demonstrates  that  the  arbitrary  multiple  a 
really  is  arbitrary,  and  is  just  as  well  chosen  as  zero.  In 
the  degenerate  case,  the  multipliers  di  and  d2  must  be  the 
same,  due  to  the  specific  relationship  between  the  repeated 
eigenvector  and  the  generalized  eigenvector  defined  in 
Eq(31).  Whereas  in  the  case  of  a  pair  of  complex  conjugate 
eigenvectors,  the  multipliers  are  chosen  to  be  the  same  for 
the  convenience  of  maintaining  conjugate  pairs  of 
eigenvectors . 

Therefore,  the  symplectic  normalization  for  each  pair 
of  eigenvectors  follows  the  same  procedure.  First  the 
matrix  of  values  F*ZF  are  found.  Then  multipliers  are  found 
by  applying  Eq(28)  co  each  pair  of  eigenvectors. 

The  value  for  a  (F’ZF)i^i+i  pair  will  commonly  be  a 
negative  or  a  complex  number.  This  inevitably  results  in  a 
set  of  symplectic  eigenvectors  that  are  complex  and  thus  a 
set  of  complex  modal  variables.  This  becomes  very 
inconvenient  for  analysis  of  the  modal  variables. 
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3.3  Real  Valued  Symplectic  Eigenvectors 

Looking  at  the  new  variational  Hamiltonian  as  defined 
in  Eq(23),  the  matrix  S  can  be  found  with  the  aid  of  Eq(22), 
since 

S  =  Z^J  =  -ZJ  (36) 


where  again,  S  is  defined  as  d^K/dY^,  or  the  second  partial 
of  the  new  Hamiltonian  with  respect  to  the  new  state  vector. 
It  will  be  the  variational  Hamiltonian,  SR,  rather  than  the 
system  Hamiltonian,  K,  that  will  be  of  concern  in  the 
following  sections.  Therefore,  the  value  of  S  will  be  of 
primary  importance. 

3.3.1  Real  Transformations  for  Non-degenerate  Cases 

For  any  system  without  a  degenerate  mode,  the  J  matrix 
will  be  of  the  form 


to^  0 
0  -co^ 


(37) 


Using  Eq(37)  in  Eq(36),  S  is  of  the  form 


0  -Oj 

co^  0 


(38) 


The  variational  Hamiltonian  of  Eq(23)  can  then  be  written  as 
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i=odd  integers 


(39) 


Sft 


=  =  5]  co,y,y 

^  1-1 


i*i  ' 


for  each  pair  of  Poincare  exponents.  Equation  (39)  provides 
the  following  equations  of  variation 


yi 


yiM 


391 

_  391 


=  “lyi 
=  -“lyi*! 


(i=odd  integers) 
{i=odd  integers) 


(40) 


which  have  the  simple  solutions 


yi  =  yi  e  ' 


yin  =  yin  e 


-a,t 


(41) 


Again,  this  is  easily  applied  when  coi  are  real,  but  for 
imaginary  <Oi  the  result  is  a  pair  of  complex  valued  modal 
variables.  This  results  in  the  need  for  another 
transformation  to  ensure  real  modal  variables.  A 
transformation  will  also  be  required  to  produce  real 
symplectic  eigenvectors  in  order  to  make  integrations  of  the 
eigenvectors  easier  to  implement  in  computer  code.  In  some 
cases,  these  transformations  can  be  accomplished  with  the 
same  transformation  matrix. 

The  matrix  that  will  transform  the  imaginary  pairs  of 
coi  in  the  J  matrix  is  defined  as 
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-  1 
'  yiL-* 


(Strang,  1988:298).  Recall  that  any  transformation  to  a 
canonical  system  must  obey  Eq(18);  therefore  in  this 
transformation,  T*ZT  must  equal  Z.  To  meet  this 
requirement,  Eq(42)  is  found  to  be 


T.  = 


1  [1  l' 


and  will  be  defined  as  a  type  A  transformation  matrix. 
Using  Tj^,  in  Eq(ll)  becomes 

*t>  =  T,E-^ 


The  matrices  E,  J,  and  S  are  also  redefined  as 


E'  =  E{TJ-^ 


0  i©^ 

-l©_,  0 


1©J  0 

s'  =  ^ 


0  I©, 
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Although  this  transformation  creates  a  real  J  and  therefore 
a  real  S,  the  resulting  symplectic  eigenvectors  will  become 
either  purely  real  or  purely  imaginary.  Ensuring  that  these 
eigenvectors  are  real  can  be  accomplished  through  another 
transformation  matrix,  or  a  type  B  transformation  defined  as 


0  I 
I  0 


(48) 


Following  the  same  method  as  with  a  type  A  transformation, 
E,  J,  and  S  become 

E"  =  E'  (Tb)  '"  (49) 


J 


n 

1 


=  t.j'at,) 


-1 


0  -Ico 
I©  0 


(50) 


S 


n 

1 


-zj'l 


-I©  0 

0  -I© 


(51) 


A  fully  populated  Tj^  matrix  will  be  block  diagonal  with  a 
for  each  pair  of  imaginary  ©^  and  an  identity  matrix  for  all 
other  pairs  of  ©i.  Likewise,  a  fully  populated  T*  matrix 
will  be  block  diagonal  with  a  for  each  pair  of  imaginary 
Pi  and  an  identity  matrix  for  all  other  pairs  of  p^. 

3.J.J  Real  Transformations  for  Degenerate  Cases 
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In  the  degenerate  case,  the  submatrix  will  always  be 
real,  so  a  type  A  transformation  will  never  be  necessary. 

As  in  the  non-degenerate  case,  the  symplectic  pi  can  be 
purely  imaginary  as  a  result  of  symplectic  normalization, 
and  therefore  a  type  B  transformation  should  be  applied. 

The  effect  on  E  will  be  the  same  as  that  in  Eq(49),  but 
since  the  original  form  of  J  and  S  in  the  degenerate  case 
are 


b  i' 

O 

,o 

= 

and  S,  = 

0  0 

1 

[o  ij 

(52) 


the  transformation  results  in 


J'l  = 


0  0 
1  0 


and  s' I 


-1  0 
0  0 


(53) 


3.4  Modal  Variables  and  E  After  Real  Canonical 
Transformation 

While  Eq(23)  is  still  valid  for  the  variational 
Hamiltonian,  the  specific  form  for  the  ith  portion  of  Eq(38) 
will  now  have  six  different  possibilities.  These  six 
different  variations  are  outlined  in  Appendix  A. 

The  final  results  of  the  last  several  sections  is  that 
Floquet  theory  can  be  applied  to  any  time  periodic  system 
and  result  in  a  real  canonical  transformation  to  modal 
variables  through  ifc(t)  =  E(t) J(t) .  There  will  also  be  a 
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Jordan-like  normal  form  J  that  is  real  valued  over  the 
period  of  the  system  (Wiesel  and  Pohlen,  1992:14). 

Regardless  of  the  specific  form  for  the  variational 
equations  of  motion,  the  time  periodic  symplectic 
eigenvectors  can  be  characterized  over  one  period  through 
Eq(14)  and  will  remain  real  over  that  period. 

In  characterizing  the  motion  of  the  system,  the  initial 
state  vector,  X,  and  the  transformation  matrix,  E,  can  be 
integrated  over  one  period  with  their  values  sampled  at 
regular  intervals.  These  values  are  then  converted  into  a 
set  of  Fourier  coefficients  to  be  used  in  construction  of 
the  perturbation  solutions  (Brower  and  Clemence,  1961:108- 
112) .  "The  advantage  of  the  Courier  representation  is  that 
the  coefficients  may  be  reassembled  into  the  periodic  orbit 
at  any  time  necessary"  (Ross,  1991:28). 

3.5  The  Restricted  Three-Body  Problem 

With  the  theory  for  canonical  Floquet  theory  developed, 
it  becomes  important  to  examine  a  specific  periodic  system 
to  validate  the  usefulness  of  the  theory.  In  the  case  of 
this  study,  the  restricted  three-body  problem  presented  by 
Ross  will  be  used.  The  particular  system  examined  is  the 
Sun-Jupiter  system.  Only  the  pertinent  equations  and 
information  will  be  presented  here,  while  the  reader  is 
referred  to  Ross  (1991)  for  detailed  derivations. 

The  definition  of  the  restricted  three-body  problem  was 
first  presented  by  Euler  in  1772.  The  problem  is  defined 
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as : 


Two  bodies  revolve  around  their  center  of  mass  in 
circular  orbits  under  the  influence  of  the  mutual 
gravitational  attraction  and  a  third  body  (attracted  by 
the  previous  two  but  not  influencing  their  motion) 
moves  in  the  plane  defined  by  the  two  revolving  bodies. 
The  restricted  problem  of  three  bodies  is  to  describe 
the  motion  of  this  third  body.  (Szebehely,  1967:8) 

In  setting  up  the  three-body  problem,  Ross  first  defines 

several  non-dimensional  quantities.  These  are: 


_ _ _  , 


M, 


«  0 


(54) 


s 


1 


5: 

' 


5a 

^7^ 


(55) 


where  and  are  the  masses  of  the  primary  bodies  and  M3 

is  the  mass  of  the  third  body  which  is  negligible.  Sj  and 
Sj  are  the  distances  from  the  primary  bodies  to  the  system 
center  of  mass.  Since  Si+S2  =  mi+m2  =  1,  these  quantities 
are  then  redefined  with  a  single  dimensionless  variable  H. 

S3  =  /Rj  =  H.  ,  S2  =  Rii  =  l-)l  (56) 

The  relationship  of  these  dimensionless  quantities  are  shown 
in  Figure  1  (Ross,  1991:4-6). 

The  Hamiltonian  for  the  restricted  three-body  system  is 
then  defined  by  Ross  as 
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(57) 


3.5.1  Periodic  Orbits  and  the  Equations  of  Variation 

The  definition  for  a  periodic  orbit  is  simply  stated  as 
the  state  vector  at  time  t  is  equal  to  the  state  vector  at 
the  initial  time,  or 

X(0)  =  7r(T)  (61) 
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Using  the  periodicity  of  a  system,  a  solution  for  the 
initial  condition  can  be  found  by  "iteratively  narrowing  the 
difference  between  the  initial  and  final  conditions"  (Ross, 
1991:9).  In  order  to  make  the  appropriate  adjustments  to 
each  iteration,  the  behavior  of  nearly  periodic  orbits  must 
be  known.  The  equations  of  variation  describe  the  motion  of 
these  nearby  trajectories.  As  defined  by  Eqs(4),  (5),  and 

(6) ,  the  equations  of  variation  for  the  three-body  problem 
are  given  as 


'o  1  1  oip^r 

-H,,  0  1  bp, 

-10  0  1  6g, 

-H,,  -1  -H,,  Oj  6p3 


A(t)X 


(62) 


These  variational  equations  of  motion  form  a  set  of  time- 
varying,  linear,  differential  equations  which  then  follow 
Eqs(7)  and  (8)  in  forming  a  fundamental  set  of  solutions,  O. 
Equations  (60)  and  (62)  can  be  numerically  integrated  to 
form  a  solution  to  Eq(7)  in  the  form  of 
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(64) 


Sgi  (X)  ■ 

^11  4^12  ^13  4*14 

6gi(0)' 

6pi  (X) 

*^21  4*22  4*23  4*24 

6p,(0) 

(X) 

4*31  4*32  4*33  4*34 

6g,  (0) 

Spj  (X) 

4*41  4*42  4*43  4*44 

6p,(0) 

L  J 

3.5.2  Determination  of  Initial  Conditions 

In  finding  a  periodic  orbit,  any  set  of  initial 
conditions  can  be  inserted  into  Eq(59)  and  integrated  for  a 
desired  x.  The  error  in  the  final  state  vector  is  the  left 
hand  side  of  Eg  (64)  and  after  inverting  the  integrated  <I>, 
the  error  in  the  initial  state  vector  is  found. 

Since  the  Hamiltonian  state  variables  are 
interdependent,  this  method  of  iteration  is  overly 
cumbersome  in  solving  for  the  initial  conditions.  The  state 
variables,  Eq(59),  can  be  shown  to  be  functions  of  H,  \l, 
qi,  and  q2  only  (Ross,  1991:17-18).  It  makes  sense  to 
maintain  the  conserved  quantity,  H,  constant  since  it  will 
allow  for  better  comparisons  in  the  perturbation  study.  The 
quantity  |A  is  system  dependent  and  therefore  constant.  This 
leaves  us  with  only  the  coordinates  qj  and  qz  to  manipulate. 
After  some  examination  it  is  noted  that  selection  of  q2=0 
also  yields  Pi=0,  or  in  other  words,  the  orbit  starting 
point  will  be  on  the  qi  axis  in  a  motion  perpendicular  to 
the  qi  axis.  The  result  of  this  selection  for  initial 
conditions  is  that  Eq(64)  can  be  reduced  to 
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(65) 


6Pi  (T) 

4>21 

^24 

6gi(0) 

(t) 

<t'34 

6p2(0) 

(Ross,  1991:13). 

But  since  the  selection  of  P2  is  not  arbitrary,  another 
revision,  as  outlined  in  discussions  with  Dr.  Wiesel, 
replaces  Pz  with  another  selectable  initial  condition,  the 
integration  time.  This  selection  then  yields  another 
quantity,  the  period,  that  will  be  as  important  as  the 
initial  coordinates  and  momenta.  Equation  (65)  can  now  be 
modified  to 


3pl (T) 

dPiCc) 

8Pi  (T) 

dt 

(0)  ■ 

6g2  (T) 

Bq2(x) 

6t 

dqJO) 

Bx 

(66) 


where 


3Pi  {X) 

Bp^(x) 

5Pi  (T)  i 

Bp^iO) 

Bq[{0) 

(0) 

dp,  (0) ; 

direct  -^2  '  '  1 

^  Bq[{0) 

3Pi('c) 


=  Pi 


(67) 


and  likewise 
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(68) 


3g2  ('') 


Pi 

<J2 


dg^  (t) 
3t 


=  P2 


The  non-direct  term  in  the  first  equations  of  Eq(68)  and 
(69)  are  due  to  the  coupling  effect  of  q^  and  Pz.  Through 
inversion  of  Eq(67)  the  required  changes  in  q^  and  the 
period  can  now  be  found- 
3.5.3  Surface  of  Section 

Finding  an  initial  guess  for  q^  and  x  is  made  easier 
through  the  use  of  surface  of  section  plots.  These  plots 
define  regions  of  periodic  motion  by  plotting  only  the 
apogee  and  perigee  points  along  a  periodic  trajectory.  A 
collection  of  these  plots  was  compiled  into  Jefferys'  An 
Atlas  of  Surface  of  Section  for  the  Restricted  Problem  of 
Three  Bodies.  As  an  example  of  a  surface  of  section  plot. 
Figure  2  shows  a  processing  elliptical  trajectory.  Given 
enough  time,  a  plot  of  only  the  perigee  and  apogee  points 
for  an  infinite  set  of  ellipses  will  form  a  set  of  rings  as 
seen  in  figure  3.  Not  all  periodic  regions  will  form  a  pair 
of  rings,  but  it  is  the  easiest  form  to  understand.  To 
examine  other  periodic  regions  the  reader  is  directed  to 
Jefferys  (1971) .  These  regions  can  then  be  searched  for 
specific  periodic  trajectories. 

Since  the  development  of  the  equations  of  motion  in 
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Jefferys  differed  from  that  of  Ross  (as  adapted  from 
Szebehely) ,  Ross  describes  the  transformation  between  the 
two  derivations  (Ross,  1991:14-17).  The  purpose  of  this 
transformation  was  so  that  the  initial  conditions  used  to 
create  Jefferys  atlas  could  be  used  as  input  for  Szebehely' s 
equations  of  motion.  The  resulting  surface  of  sections 
could  be  directly  compared  to  the  atlas  developed  by 
Jefferys . 

For  the  Sun- Jupiter  system,  11=0.00095388.  With  the 
selection  of  3.15  for  the  Jacobian  constant,  J,  the 
resulting  surface  of  section  can  be  seen  in  figure  4.  The 
relationship  between  the  Jacobian  constant  and  the  system 
Hamiltonian  is 

if  =  -J]  (69) 


(Ross,  1991:16). 

3.6  Perturbation  Theory  on  the  Restricted  Three-Body 
Problem 

With  the  assistance  of  the  surface  of  section  plot  and 
the  iterative  integration  of  the  equations  of  motion,  a 
periodic  trajectory  is  found.  The  next  step  is  to  find  out 
what  effect  perturbations  to  the  initial  conditions  have  on 
the  stability  of  the  orbit. 

In  examination  of  the  effect  of  perturbations  to  the 
initial  conditions,  two  different  representations  of  the 
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Figure  2;  An  Elliptical  Trajectory  Precessing 
About  Primary  l-ji  (Ross,  1991:19) 
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Figure  3:  The  Surface  of  Section  of  the 
Elliptical  Trajectory  (Ross,  1991:20) 
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cross 


J  =  3.15 


=  0.00095388 


Figure  4:  The  Surface  of  Section  for  the 
Sun-Jupiter  System 
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nearly  periodic  orbit  will  be  constructed.  Both 
representations  will  use  the  <I>  found  from  the  integration  of 
the  periodic  initial  conditions  for  one  period.  The  <I> 
matrix  will  be  transformed  into  a  real  valued,  symplectic 
eigenvector  matrix,  E,  and  a  Jordan  normal-like  matrix,  J, 
of  Poincard  exponents. 

3.6.1  The  Exact  Representation 

In  the  first  representation,  the  periodic  X  and  E  will 
be  integrated  over  one  period  and  an  evenly  spaced  sampling 
of  these  quantities  will  be  converted  into  a  set  of  one 
hundred  Fourier  coefficients.  A  nearly  periodic  trajectory, 
close  to  the  original  trajectory,  will  then  be  integrated 
over  time.  The  integration  time  required  for  the  nearly 
periodic  trajectory  will  be  many  times  the  orbital  period  of 
the  system  and  is  on  the  order  of  the  period  of  an 
oscillatory  Poincard  exponent.  Subtracting  the  periodic 
trajectory  from  the  nearly  periodic  defines  an  x(t)  (as 
described  in  Eq(6)),  and  using  the  canonical  Floquet 
transformation, 

y(t)  =  E-Mt)X(t)  (70) 

creates  a  set  of  modal  variables  over  the  integration  time. 
This  set  of  modal  variables  can  then  be  plotted  as  an  exact 
representation  of  the  nearly  periodic  orbit  (Ross,  1991:24). 

3.6.2  The  Expanded  Representation 

For  the  second  representation,  the  original  Hamiltonian 
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is  canonically  transformed  into  the  modal  variables  and 
expanded  in  a  Taylor  series  which  produces 


d^H(y) 


(71) 


YiYj^ 


ly*o 


where  J=0  centers  the  expansion  on  the  periodic  trajectory. 
In  tensor  notation  Eq(71)  becomes 

(72) 

K(y)  =  HCU)  . 


The  first  term  in  the  expansion  is  the  Hamiltonian 
for  a  periodic  orbit,  and  is  a  constant.  The  second, 
or  linear  term  is  identically  zero,  because  it 
describes  the  motion  of  the  periodic  trajectory  with 
respect  to  itself.  The  third,  or  quadratic  term  is  the 
Floquet  problem,  and  becomes  a  constant  coefficient, 
linear  system  in  the  new  variables.  Since  the 
magnitude  of  the  modal  state  vector  is  small  compared 
to  one,  the  expansion  is  truncated  after  the  fourth 
term.  (Ross,  1991:31-32) 

The  new  modal  Hamiltonian  then  becomes 


fC(p)  =  H{XJ)^^y^Sy*^H,^,[V)y,yjy,  (73) 


where  the  ith  portion  of  the  quadratic  term  will  follow  one 
of  the  forms  outlined  in  Appendix  A. 

The  third  order  tensor  can  be  expanded  to  the  very 
cumbersome  form  of 
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(74) 


*  c.yiVa  ^  C3yiy3  ^  c.y^y^  ^  c^y^yl 

+  ^.yiyaya  ^  c^y.y^y,  +  c^y^Ys  +  c^yiy3y4  *  c^oy^y' 

+  c^.yi  +  c,,yiy^  +  c.^y^y,  +  c^^y^ys"  +  c^sy^yjy, 

+  Ci^y^y,"  +  c.^yj  +  c^gYaY,  +  c,,y3y4  +  c^oY,^ 

where  the  coefficients  c^  are  defined  by  the  triple 
sununation  i=l  to  4,  j=l  to  4,  and  k=l  to  4  over  the  right 
hand  sides  of 
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Cl  (t) 

= 

C,(t) 

= 

llP JlP k2*P llP J2P kl'*'P I2P ]lP kl^ 

C^(t) 

= 

-^^l]k  [PllPjlP*3"''PllPj3P*l"'’Pl3P:;2Pjti] 

C4  (t) 

= 

^P llP JlP k4'^P llP jlP kl'^P  14P jlP kl^ 

Cs  (t) 

= 

-^^l]k  ^P llP J2P k2'’'P jlP k2'*'P  12P j2p kl^ 

C6(t) 

= 

Jk  t  PilP J^P^-^PuP J3p*2+Pi2p ]lPk3 

'*'P  12P  :}3pkl'^P  13P]lPk2'*'P  13P  J2Pkl 

C7  (t) 

= 

-^f^l]k^PllP]2Pk4^PllP]4Pk2^Pl2p]lPk4 

"^P 12P  J4P  kl'*'P  14P  JlP  k2'*'P  14P  ]2P  kl 

C8(t) 

= 

^ijk  [ P iiP  jsP )c3^P  13P  jiP ;t3''’P  13P :33P Jti  1 

C,{t) 

= 

-^^l]k^P  ilP  ]3Pk4'*'P  llPj4Pk3'*'P  13P  ]lPk4 

'*'P  13P  J4Pkl'*'Pl4P  ]lPk3'*'P  14P  J3Pkl 

Cio(t) 

= 

^P  llP  j4P  k4'^P  14P  JlP  k4'*’P  I4P  ]4P  kl^ 

Cii(t) 

s 

-^^IjkP  12P  ]2Pk2 

Ci2(t) 

s 

~^^l]k^P  I2P  J2P  k3*P  I2P  j3p  Jt^'^’P  J3P  J2P  k2  ^ 

Ci3(t) 

= 

-g ^Ijk  t  P 12P  j:Pk4 '*'P  12P  J4P )C2''‘P  14P  32P k2 ^ 

Ci4  (t) 

= 

-g^lj*  ^P j:P k3'*'P  )2P k3'*'P  13P j3p k2^ 

Ci5(t) 

= 

-^^l]k  ^P 12P ]3P k4'*‘P I2P ]4P ks'^P I3P J2P k4 

'*'P  13P  j4P  k2'*'P  14P  ]2P  k3'*'P  14P  ]3P  k2 

Ci6(t) 

= 

-^Hl]k^Pl2P}4Pk4^Pl4Pj2Pk4*Pl4P:}4Pk2'i 

Ci7  (t) 

= 

-g  JkP  13P  J3P  fc3 

Cie(t) 

= 

-g^lj*  ^P 13P  j3P k4'*'P  13P J4P k3^P  14P JsP k3^ 

Ci9(t) 

= 

-^^Ijk  ^P I3P j4P k4'*'P  14P ]3P k4'*'P I4P j4P k3^ 

C2o(t) 

= 

^^l]kPl4p]4Pk4 

Here,  the  variable  p^,  is  the  element  in  the  /nth  row  and  the 
nth  column  of  the  transformation  matrix,  E.  The  periodic  St 
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and  E  will  again  be  integrated  over  one  period,  but  for  this 
representation,  at  each  of  the  evenly  spaced  samplings,  the 
expanded  modal  coefficients,  c^,  will  be  calculated.  Like 
the  sampled  sets  in  the  exact  representation,  each  of  the 
sampled  sets  of  c^  will  be  converted  into  a  set  of  one 
hundred  Fourier  coefficients. 

Starting  with  the  initial  modal  variables,  found  in  the 
exact  solution  at  t=0,  the  expanded  Hamiltonian  is 
integrated  over  the  same  time  as  the  exact  solution.  The 
results  can  be  plotted  as  the  expanded  solution  and  can  be 
compared  to  the  exact  solution.  The  modal  equations  of 
motion  come  from  the  expanded  Hamiltonian  and  are 


«  =  =  9<^£)  .c 


C2yi  ^2c5yiy3+C5yiy3  +  c,yjy, 


+ 3  + 2  Ci2y,y3 + 2  + c,,yy  +  c,  3y3y,  +  ,y! 
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where  the  term  d(K2)/dy^  depends  on  the  form  of  S  in  the 
quadratic  term  of  the  modal  Hamiltonian.  Appendix  A  also 
defines  the  forms  of  d(K2)/9yi  for  the  various  types  of 
eigenvector  pairs. 
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IV.  Software 


The  software  programs  used  in  this  study  are  coded  in 
Fortran  77.  All  seven  main  programs  and  eleven  subroutines 
can  be  found  in  Appendix  B.  The  majority  of  the  programs 
were  written  with  the  restricted  three-body  problem  in  mind, 
but  any  of  the  software  can  easily  be  modified  to  any 
periodic  Hamiltonian  system. 

4.1  Surface  of  Section 

The  first  program  used  in  the  study  is  the  program 
SECTION.  The  core  of  the  program  was  written  by  Ross  (1991) 
with  the  purpose  of  duplicating  the  surface  of  section  plots 
created  by  Jefferys  (1971) .  The  program  is  written  to 
integrate  several  initial  conditions,  while  storing  the 
apoapsis  and  periapsis  crossings  in  a  plot  file.  While  the 
program  is  simple  in  its  concept,  finding  initial  conditions 
that  produce  periodic  regions  in  the  phase  space  is 
guesswork  at  best.  Validation  of  this  program  is 
accomplished  through  direct  comparison  with  the  atlas  of 
surface  of  sections  created  by  Jefferys. 

4.2  Determination  of  Periodic  Initial  Conditions 

Once  the  surface  of  section  plot  has  been  created, 

determination  of  exact  initial  conditions  for  a  periodic 
trajectory  must  be  found.  The  reader  should  recall  that 
this  was  done  by  selecting  an  initial  condition  where  the 
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trajectory  was  crossing  the  axis  at  the  initial  time. 
Therefore  the  initial  condition  on  qj  is  zero;  the  guessed 
initial  condition  on  q^  is  found  by  estimating  the  center  of 
a  periodic  region  on  the  qi  axis.  The  third  initial 
condition,  the  period,  is  estimated  by  examining  the  time  it 
takes  for  a  nearby  orbit  to  nearly  return  to  its  initial 
location.  These  estimated  initial  conditions  are  loaded 
into  the  program  PERIOD  which  iteratively  integrates  and 
corrects  qi  and  x  until  at  the  end  of  one  period,  the 
trajectory  has  returned  to  the  initial  qi  and  q^  position. 

A  benefit  of  this  program  is  that  it  is  self  checking,  it 
accomplishes  this  by  forcing  the  final  conditions  and 
initial  conditions  to  match.  Verification  of  the  coded 
equations  of  motion  is  the  only  check  needed  for  this 
program,  and  this  is  best  done  by  hand. 

Once  the  initial  conditions  are  found,  the  program 
extracts  the  eigenvalues,  eigenvectors,  and  Poincar6 
exponents  of  the  system. 

4.3  Symplectic  Normalization 

The  type  of  symplectic  normalization  required  is 
determined  by  examination  of  the  types  of  eigenvector  and 
Poincar6  exponent  pairs  found  in  program  PERIOD.  Some 
rearrangement  of  the  eigenvectors  and  eigenvalues  may  be 
required  in  order  to  ensure  positive/negative  or  degenerate 
pairs  are  kept  together.  Once  properly  arranged,  the 
eigenvectors  are  fed  into  program  EXSYRL  (EXTENDED, 
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SYMPLECTIC,  REAL) .  The  purpose  of  this  program  is  to  find 
an  extended  eigenvector  if  needed,  convert  the  standard 
eigenvectors  into  symplectic  eigenvectors,  and  ensure  that 
the  symplectic  eigenvectors  and  the  J  matrix  are  real 
valued.  The  best  check  of  this  program  is  to  ensure  the 
final  forms  of  E  and  J  satisfy  Eg (11)  at  several  different 
times  over  the  period  of  motion. 

4.4  Storage  of  the  Periodic  Trajectory  and  Hamiltonian 
Coefficients 

The  final  E  and  J  matrices  are  loaded  into  programs 
FLOQUET  and  HAMILTONIAN  to  create  the  sets  of  Fourier 
coefficients  for  either  the  exact  or  the  expanded 
perturbation  problems. 

FLOQUET  integrates  the  periodic  orbit  and  extracts  two 
hundred  evenly  spaced  values  of  the  state  vector  and  the 
symplectic  eigenvector  matrix.  These  sets  are  then  fed  to 
subroutine  FOURIER  which  converts  them  into  one  hundred 
pairs  of  Fourier  coefficients  allowing  the  periodic 
trajectory  to  be  reformed  at  any  time. 

HAMILTONIAN  also  integrates  the  periodic  trajectory  and 
extracts  the  state  vector  and  the  E(t)  matrix  values  at  two 
hundred  evenly  spaced  points  on  the  trajectory.  It  then 
performs  a  Taylor's  series  expansion  of  the  Hamiltonian  to 
find  the  twenty  Hamiltonian  coefficients,  c^it) ,  at  each  of 
the  two  hundred  points.  These  sets  of  Ci(t)  are  also  loaded 
into  FOURIER  to  be  converted  into  one  hundred  pairs  of 
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Fourier  coefficients  for  each  set  of  Ci(t). 

One  check  for  these  two  programs  is  to  ensure  that  the 
state  vector  and  E  have  returned  to  their  initial  condition 
after  one  period.  This  should  not  be  a  problem  since  this 
same  check  is  performed  in  PERIOD.  A  valuable  check  of  the 
Fourier  coefficients  is  performed  by  summing  the  Fourier 
cosines.  Each  cosine  sum  should  equal  the  initial  condition 
of  the  function  the  Fourier  coefficients  describe  (Brower 
and  Clemence,  1961:110). 

4.5  The  Perturbation  Solutions 

The  exact  and  the  expanded  perturbation  problems  are 
analyzed  by  programs  of  the  same  names.  Program  EXACT  takes 
the  Fourier  coefficients  created  in  FLOQUET  and  reassembles 
the  periodic  trajectory  conditions  on  demand.  Perturbed 
initial  conditions  are  integrated  in  time,  and  at  regular 
intervals  the  state  vector  is  extracted.  The  periodic  state 
is  subtracted  from  the  perturbed  state  and  the  result  is 
converted  into  the  modal  variables  by  the  corresponding 
E{t)  . 

Program  EXPANDED  requires  not  only  the  Fourier 
coefficients  created  in  program  HAMILTONIAN,  but  also  the 
initial  modal  conditions  as  found  by  program  EXACT. 

EXPANDED  then  integrates  the  modal  variables  directly.  Both 
programs  extract  the  modal  variables  at  regular  intervals 
and  send  the  results  to  plot  files.  If  the  two  programs  are 
working  correctly,  there  should  be  a  certain  amount  of 
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correspondence  of  the  plots  for  the  two  solutions.  The 
actual  amount  of  correlation  is  a  topic  for  the  next 
chapter . 
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V.  Results  and  Discussions 


5.1  Canonical  Floquet  Theory 

The  most  significant  result  of  the  canonical  Floquet 
theory  is  that  it  works.  The  application  of  canonical 
Floquet  theory  is  currently  limited  to  Hamiltonian  systems. 
There  appears  to  be  no  limits  to  the  degrees  of  freedom  that 
can  involved  as  long  as  the  system  can  be  written  as  a 
periodic  Hamiltonian.  This  now  opens  the  doors  to  a  great 
number  of  possibilities  for  perturbation  solutions  to 
periodic  systems.  Thus,  the  results  of  the  three-body 
perturbation  solution  will  be  of  more  interest  than  any 
further  discussion  of  the  theory. 

5.2  The  Three-Body  Perturbation  Problem 

Recall,  the  particular  three-body  system  of  interest  in 
this  study  is  the  Sun-Jupiter  system;  this  means  a 
|1  =  0.00095388.  The  surface  of  section  for  this  system  was 
presented  in  Figure  4.  The  first  task  was  to  find  a 
periodic  region  in  this  phase  space.  The  surface  of  section 
was  created  using  a  Jacobian  constant  of  J  =  3.15  and  a 
dozen  initial  state  conditions  along  the  x  axis.  In 
examining  the  x  axis  between  the  two  primary  masses  (i.e. 
between  -1  and  0) ,  there  appears  to  be  two  main  periodic 
regions;  one  centered  at  approximately  x  =  -0.3  and  the 
other  at  approximately  x  =  -0.5  on  the  surface  of  section. 
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The  five  crescent  shaped  structures  in  between  the  first  and 
second  central  rings  define  the  periodic  region  associated 
with  X  =  -0.3,  while  there  are  four  outer  crescents  that 
define  the  periodic  region  associated  with  x  =  -0.5.  Figure 
5  shows  a  close-up  view  of  the  x  =  -0.5  region  and  provides 
an  initial  guess  of  x  =  -0.51  for  the  periodic  initial 
condition,  the  center  of  the  crescent  structure.  After 
iteration  of  the  guessed  initial  conditions  with  program 
PERIOD,  the  actual  initial  conditions  are  determined  to  be 
X  =  -0.5136620321,  y  =  0.0,  and  a  period  of  6.256580411  (all 
rounded  to  ten  significant  figures)  for  the  given  Jacobian 
and  |1.  Converting  these  conditions  into  the  initial 
conditions  for  Szebehely's  equations  yields 
q,  =  -0.5127081521,  qj  =  0.0,  and  H  =  -1.574523515,  with  the 
same  period  of  motion.  This  particular  periodic  trajectory 
produced  a  positive/negative  imaginary  pair  of  Poincare 
exponents  as  well  as  a  degenerate  pair,  making  it  a  good 
test  case  for  the  canonical  Floquet  theory  as  well  as  the 
perturbation  analysis. 

5.2.1  Unperturbed  System 

The  starting  point  for  analysis  in  the  study  was  to 
examine  the  time  history  of  the  modal  vectors  for  the 
unperturbed  case.  Since  an  imaginary  pair  of  (o^  produces  a 
sine  and  a  cosine  function  in  the  modal  variables,  a  plot  of 
one  versus  the  other  should  produce  a  circle.  For  the 
unperturbed  case,  the  modal  vectors  should  equal  zero  for 
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Figure  5:  Close-up  View  of  Periodic  Region  of  the 
Sun-Jupiter  Surface  of  Section 
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all  time.  Figure  6  shows  the  first  and  second  modal  vectors 
plotted  against  each  other.  The  other  two  modal  vectors  are 
simply  plotted  against  time  in  Figures  7  and  8.  In  all 
three  figures,  the  scale  is  such  that  maximum  deviation  will 
be  seen,  but  in  all  cases  the  magnitude  of  the  deviation  is 
zero  for  purposes  of  this  analysis.  Looking  at  the 
deviations  in  this  detail  do  show  trends  in  the  data  that 
can  be  compared  to  data  in  the  perturbation  cases.  The  fact 
that  any  error  shows  at  all  is  almost  certainly  due  to 
truncation  error  in  all  the  calculations  made  on  this  data. 
5.2.2  Perturbed  Systems 

For  this  study  only  two  initial  values  were  perturbed 
for  examination,  first  qj  and  then  J;  J  can  be  translated 
into  a  change  in  H.  Each  quantity  affected  the  value  of  Pz 
through  the  equations  of  motion,  but  no  other  initial 
conditions  are  changed.  Physically,  a  change  in  q^,  and  Pz, 
means  a  change  in  position  on  the  x  axis  of  the  system  in 
Figure  1  along  with  a  change  in  the  y  velocity,  Pz,  in  order 
to  stay  on  the  same  Hamiltonian  surface,  H.  The  change  in 
the  Hamiltonian  value  allows  for  the  position  to  remain  the 
same,  but  a  change  in  the  y  velocity  is  still  required. 

5. 2. 2.1  Changes  in 

As  in  any  perturbation  problem,  the  magnitude  of  the 
change  must  be  small  in  order  to  keep  the  truncated, 
expanded  Hamiltonian  valid.  Since  the  distance  between  the 
two  primary  masses  has  been  non-dimensional i zed  to  a  value 
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Figure  8:  Unperturbed  Time  History  of 
Modal  Variable 


of  one,  it  would  be  accurate  to  say  that  a  perturbation  of 
0.001  for  the  third  body  in  the  system,  or  a  tenth  of  a 
percent,  is  significant.  In  physical  dimensions  this  would 
be  an  approximate  position  change  of  775,000  km,  or  more 
than  twice  the  distance  to  the  moon  from  earth.  The  initial 
change  to  the  position  was  chosen  to  be  =  le-8,  or  about 
7.75  km.  It  can  be  seen  in  Figure  9  that  the  perturbed 
solution  of  the  oscillatory  modal  variables  describes  a 
visually  perfect  circle.  Some  particular  items  of  interest: 
1)  the  circle  does  close  on  itself  as  can  be  seen  by  an 
overlapping  region  near  the  left  side  of  the  plot,  2)  both 
the  exact  and  the  expanded  solution  match  precisely  for  this 
initial  displacement,  3)  the  amplitude  of  the  and  y2 
modal  variables  is  about  three  times  the  initial  position 
displacement,  and  4)  a  time  history  of  the  first  two  modal 
variables.  Figures  10  and  11,  shows  that  these  two  variables 
do  in  fact  describe  smooth  sine  and  cosine  functions  at  this 
level  of  displacement.  The  time  history  for  the  perturbation 
in  modals  ya  and  y,  is  not  presented  at  this  point  because 
the  magnitudes  of  the  perturbations  is  still  nearly  zero. 

The  exact  and  expanded  solutions  of  y^  and  y^  match 
perfectly  at  =  le-8. 

The  magnitude  of  the  perturbation  is  then  increased  to 
6qi  =  le-6,  or  about  775  km  (the  distance  from  San  Diego  to 
San  Francisco) .  Although  this  is  an  enormous  distance  for 
an  Earth  satellite  to  be  out  of  orbit,  is  it  a  significant 
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perturbation  for  a  periodic  trajectory  in  the  Sun-Jupiter 
system?  Figure  12  shows  that  although  there  is  a  slight 
waver  in  the  circular  path  of  the  first  two  models,  it  is 
still  very  close  to  cyclic  behavior.  The  time  histories  of 
the  first  two  models  show  no  difference  to  the  sine  and 
cosine  functions  of  Figures  10  and  11.  Note  also  that  the 
amplitude  of  these  models  is  still  approximately  three  times 
that  of  the  initial  displacements  and  the  exact  and  expanded 
solutions  still  follow  each  other  perfectly.  At  this  point, 
the  magnitude  of  the  third  and  fourth  models  (Figures  13  and 
14)  are  becoming  more  significant,  but  they  are  still  at 
least  one  hundred  times  smaller  than  the  first  two  modals. 

It  is  interesting  to  note  how  closely  the  exact  and  expanded 
solutions  match  in  all  modals  at  this  point. 

Increasing  the  initial  displacement  even  further,  to 
6qi  =  5e-5,  or  about  39,000  km  (three  times  the  diameter  of 
the  earth) ,  a  drastic  change  is  found  in  not  only  the  modal 
response,  but  also  the  correlation  of  the  exact  and  expanded 
solutions.  As  can  be  seen  in  Figure  15,  although  the  exact 
solution  does  close  on  itself,  it  follows  a  very  erratic 
path  that  quickly  loses  its  circular  appearance.  The 
jagged,  sharp  changes  in  many  of  the  figures  at  this  point 
are  partially  do  to  the  large  steps  between  points  in  the 
plot,  but  most  of  the  problem  is  occurring  because  the 
periodic  trajectory  is  dissolving.  This  can  be  seen  in 
Figure  16  where  the  plot  of  the  oscillatory  modals  in  the 
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Figure  12:  Oscillatory  Modal  Variables, 
Yi  versus  y2,  Sqj  =  le-6 
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Figure  15;  Oscillatory  Modal  Variables, 
Yi  versus  yj,  Sq^  =  5e-5 
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expanded  solution  is  shown  alone.  The  integration  time  is 
no  longer  enough  to  allow  these  modal  variables  to  close  on 
themselves.  The  time  histories  of  the  four  modal  variables 
for  both  the  exact  and  the  expanded  solutions  are  shown  in 
Figures  17  to  21.  In  every  case,  the  expanded  solution 
appears  to  be  lagging  the  exact  solution.  A  larger 
deviation  in  form  also  starts  showing  up  in  modals  Yj  and  y4 
of  Figures  20  and  21.  It  appears  obvious  that  even  at  this 
point,  the  expanded  solution  is  no  longer  accurate.  The 
periodic  nature  of  the  trajectory  does  seem  to  hold  for  the 
exact  solution,  but  with  such  erratic  behavior  over  part  of 
the  time  history  it  is  doubtful  that  the  third-body  in  this 
system  could  maintain  its  trajectory  for  long. 

Finally,  as  the  change  to  continues  to  increase,  the 
expanded  solution  loses  all  its  coherence  first,  as  can  be 
seen  in  Figure  22,  while  the  exact  solution  becomes  entirely 
distorted  by  =  le-3  in  Figure  23.  Looking  more  closely 
at  the  first  two  modal  variables  when  6qi  =  le-3  (twice  the 
distance  to  the  moon) ,  it  can  be  seen  that  the  y^  modal  has 
completely  lost  the  cosine  shape  and  looks  a  great  deal  like 
the  y3  modal  inverted.  The  Y2  modal  has  lost  the  smooth 
sine  shape  and  appears  to  be  approaching  the  same  shape  as 
the  y^  modal  (see  Figures  24  and  25) .  The  shape  of  the  Ya 
and  y4  modals  has  not  changed  since  the  initial 
perturbation,  while  the  magnitude  has  increased  with  the 
perturbations.  Meanwhile,  the  expanded  solution  for  the 
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Figure  18:  Time  History  of  Modal 
Variable  Expanded  Solution,  Sqj  =  5e-5 
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Figure  20:  Time  History  of  Modal 
Variable  yj,  6qi  =  5e-5 
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Figure  21:  Time  History  of  Modal 
Variable  y^,  6qi  =  5e-5 


67 


modal 


J  =  3.15 

=  0.00095388 


6J 

6qi 


0.0 

1 .  Oe-4 


-  Exact 

XXX  Expanded 


0 


CNJ 


0.0005 


-0.001 


0 


0.0005 


modal  yi 


Figure  22:  Oscillatory  Modal  Variables, 
Yi  versus  yj,  6qi  =  le-4 
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Figure  23:  Oscillatory  Modal  Variables, 
Yi  versus  yj,  Exact  Solution,  6qi  =  le-3 
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and  modals  has  begun  to  lose  coherence  and  the  value  of 
these  modals  is  above  the  double  precision  ability  of  thf^ 
computer  in  both  cases.  In  either  solution,  the  system  is 
no  longer  periodic. 

5. 2. 2. 2  Changes  in  J 

As  with  the  6qi  cases,  the  6J  cases  will  start  with 
6j  =  le-8.  The  change  in  J,  or  H,  is  not  as  easy  to 
conceptualize  as  the  change  in  q^,  but  the  results  are  just 
as  interesting.  Figure  26  shows  that  the  and  y.^  modals 
again  form  a  visually  perfect  circle  at  the  initial  change 
in  the  Jacobian/Hamiltonian  and  again,  the  exact  and 
expanded  solution  match  perfectly.  The  first  two  modals 
still  describe  their  sine  and  cosine  waves  and  are  at  this 
point  still  smooth  in  structure.  Meanwhile,  the  modals  y^ 
and  y4  are  still  too  small  to  draw  any  significant 
conclusions.  As  the  change  in  the  Jacobian  value  is 
increased  to  6J  =  le-6,  a  similar  distortion  to  the  circular 
plot  of  yi  and  yj  begins  to  appear.  This  distortion  can  be 
seen  in  Figure  27  for  the  two  modals  plotted  versus  each 
other  and  in  Figures  28  and  29  for  the  time  history  of  the 
modals.  Modals  y3  and  y^  are  also  shown  in  Figures  30  and 
31  for  the  same  change  in  the  Jacobian.  At  this  point, 
there  are  two  notable  differences  in  the  modals  with  a 
change  in  the  Jacobian  as  compared  to  those  with  a  change  in 
qi .  First,  modals  y^  and  y^  tend  to  be  more  erratic  in  the 
earlier  part  of  their  time  history.  It  is  surmised  that 
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Figure  26;  Oscillatory  Modal  Variables 
Yi  versus  yj,  6J  =  le-8 
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Figure  27:  Oscillator^'  Modal  Variables 
Yi  versus  Yj,  6j  =  le-6 
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Figure  29:  Time  History  of  Modal 
Variable  y2/  5J  =  le-6 
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this  increased  erratic  behavior  at  the  beginning  is  damped 
out  later  due  to  the  influence  of  the  primary  masses.  The 
result  is  a  change  in  the  trajectory  from  a  less  stable  to  a 
more  stable  state  for  the  new  value  of  the 
Jacobian/Hamiltonian.  Secondly,  the  y4  modal  continually 
increases  in  magnitude  throughout  the  integration  period;  in 
fact,  it  is  exactly  a  linear  increase  to  the  modal  value. 

The  most  likely  cause  of  this  is  that  this  modal  condition 
is  associated  with  the  time  along  the  trajectory,  and  the 
change  in  the  Hamiltonian  value  causes  a  constant  drift  in 
the  time  location. 

The  change  in  the  Jacobian  is  then  increased  to 
6j  =  le-5  and  5J  =  5e-5,  and  several  new  deviations  are 
observed.  In  Figures  32  and  33,  the  modals  y^  versus  y2, 
are  notably  not  converging  on  each  other  at  the  initial 
condition  (located  at  nine  o'clock  on  the  plot),  unlike  the 
case  in  all  the  y^  versus  y2  plots  for  the  Sqj  cases.  Yet, 
like  the  6qi  cases,  there  are  still  two  "calm"  regions  in 
the  path  of  the  modal  circles.  The  first  occurs  just  above 
the  three  o'clock  point  on  the  path  and  the  second  occurs  at 
about  ten  o'clock.  In  comparison,  looking  at  Figures  15  and 
22,  even  in  the  highly  perturbed  cases  these  regions  appear 
at  exactly  three  and  nine  o'clock.  A  similar  pattern  is 
seen  when  comparing  the  yj  n'.'^dal  time  histories  in  Figures 
34  and  35  to  those  of  Figures  17  and  24.  This  would  also 
support  the  idea  that  there  is  a  drift  along  the  time  domain 
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Figure  33;  Oscillatory  Modal  Variables, 
Yi  versus  Yir  =  5e-5 
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Figure  35:  Time  History  of  Modal 
Variable  yi.  Exact  Solution,  6J  =  5e-5 


of  the  system  and  that  it  is  affecting  all  the  modal 
variables.  Also  note,  that  through  6j  =  le-5,  the  expanded 
solution  has  closely  followed  the  exact  solution.  But,  at 
6j  =  5e-5,  the  expanded  solution  again  falls  off. 

Finally,  the  change  in  the  Jacobian/Hamiltonian  is 
increased  to  6J  =  le-3,  where  the  same  trends  noted  before 
are  amplified  in  magnitude.  A  close  look  at  Figures  36 
through  40  shows  that  the  second  erratic  region  for  each 
modal  variable  is  shorter  and  smaller  in  magnitude  than  the 
first,  and  it  appears  that  a  third  region  is  starting.  If 
these  regions  continue  to  get  shorter  and  smaller,  it  would 
support  the  theory  that  the  perturbed  trajectory  is  settling 
into  a  new  trajectory  that  may  be  more  stable  than  the 
original . 

5 . 3  Summary 

In  the  end  analysis,  it  appears  that  the  exact  solution 
tracks  the  modal  variables  slightly  better  at  larger 
perturbations  than  the  expanded  solution.  Since  the 
correlation  of  the  two  solutions  is  so  accurate  at  small 
perturbations,  this  suggests  that  the  biggest  fault  in  the 
expanded  solution  is  the  truncation  of  the  Hamiltonian  after 
the  cubic  term.  The  nearly  perfect  comparison  of  the  two 
solutions  also  lends  credibility  to  both  methods  used  in 
this  study  for  the  perturbation  problem.  Because  the 
perturbation  problem  was  secondary  to  this  study,  the 
interpretation  of  the  perturbed  solutions  presented  may  not 
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Figure  36:  Oscillatory  Modal  Variables, 
Yi  versus  yj,  Exact  Solution,  6J  =  le-3 
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Figure  40:  Time  History  of  Modal 
Variable  Exact  Solution,  6J  =  l€ 


be  the  only  explanation  and  the  reader  is  encouraged  to  make 
their  own  examination  of  this  and  other  periodic  solutions 
using  the  canonical  Floquet  theory. 
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VI.  Conclusions  and  Recommendations 


The  major  conclusion  of  this  study  is  that  Floquet 
theory  can  be  used  for  canonical  transformations  to 
real-valued  modal  variables  for  use  in  perturbation  research 
on  periodic  systems.  Any  transformation,  T,  to  the 
eigenvectors  or  exponents  of  the  Floquet  solution  must  obey 
the  relation  T*ZT  =  Z  in  order  for  the  transformation  to  be 
a  canonical  one.  Special  considerations  must  be  made  for 
each  Poincard  exponent  pair  type,  with  extra  emphasis  on  the 
degenerate  case  of  zero  Poincard  exponents. 

In  support  of  the  major  finding,  the  restricted  three- 
body  solution  was  analyzed  with  both  an  exact  and  an 
expanded  solution  in  the  modal  variables.  The  results 
showed  that  the  exact  and  expanded  solutions  agreed 
extremely  well  for  very  small  perturbation.  At  larger 
perturbations,  as  expected,  the  expanded  solution  was  no 
longer  accurate  due  to  the  truncation  of  the  expansion  after 
the  cubic  term.  The  exact  solution,  while  remaining 
graphically  periodic  at  higher  perturbations,  also  exhibited 
irregular  displacements  at  these  higher  perturbations  and 
would  likely  result  in  complete  lose  of  the  periodic 
reference  trajectory  over  time. 

The  recommended  follow-on  to  this  work  would  involve  a 
much  more  detailed  analysis  of  the  three-body  problem  to 
include  expanding  the  analysis  to  a  three-dimensional 
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problem.  Mass  systems  other  than  the  Sun- Jupiter  system 
would  create  more  interesting  surface  of  section  plots  and 
would  doubtlessly  result  in  much  more  interesting  modal 
analyses.  There  is  also  the  need  for  using  canonical 
Floquet  theory  on  other  periodic  systems,  not  only  to 
further  validate  this  approach  to  canonical  transformations, 
but  to  assess  the  full  usefulness  of  this  approach  in 
perturbation  problems.  Another  system  that  would  be  of 
particular  interest,  because  of  its  requirement  for  precise 
operation,  would  be  that  of  rotating  blades  on  helicopters 
and  jet  engines. 
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Appendix  A:  Six  Modal  Equations  of  Variation  Types 


All  six  cases  for  the  modal  equations  of  variation  are 
listed  in  the  following  pages.  Note  that  the  df3i/dy^  is 
equivalent  to  3(K2)/3yi  term  for  the  expansion  of  the 
variational  modal  Hamiltonian. 
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Case  1:  Poincar^  exponents  are  a  positive/negative  pair 
of  reals,  and  the  symplectic  eigenvectors  are  real.  The 
initial  and  final  forms  of  J  are 


CO  0 
0  -CO 


(A-1) 


the  resulting  S  is  then 


(A-2) 


and  the  modal  Hamiltonian  SR  is 

^  (yiy^co+y.y^co)  =  y^y^co  (A-3) 


Finally,  the  equations  of  variation  are 


Yi 

y2 


d9t 

_  asR 


=  y,(0 
-  -y^co 


(A-4) 
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Case  2:  Poincar^  exponents  are  a  positive/negative  pair 
of  reals,  and  the  symplectic  eigenvectors  are  imaginary. 

The  forms  of  J  are 


(A-5) 


and  the  modal  Hamiltonian  9?  is 

91  =  ^y*S"y  =  =  -yiyzto  (A~7) 


Finally,  the  equations  of  variation  are 


yi  = 
72  = 


=  -yi© 

=  y^© 


(A-8) 
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Case  3:  Poincare  exponents  are  a  positive/negative  pair 


of  imaginaries,  and  the  symplectic  eigenvectors  are 
The  forms  of  J  are 


col  0 
0  -col 


J' 


0  -CO 
CO  0 


and  the  resulting  S'  is 


S'  =  -ZJ' 


-CO  0 
0  -CO 


and  the  modal  Hamiltonian  91  is 

9f  =  ^y^S'y  =  --^co(y^+y2) 


Finally,  the  equations  of  variation  are 


dSt 

Vi  = 

dSt 


real . 


(A-9) 


(A-10) 


(A-11) 


(A-12) 
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Case  4:  Poincar^  exponents  are  a  positive/negative  pair 
of  imaginaries,  and  the  symplectic  eigenvectors  are 
imaginary.  The  forms  of  J  are 


col  0 
0  -col 


0  -CO 
00  0 


J" 


0  CO 
-CO  0 


(A-13) 


and  the  resulting  S' '  is 


S" 


-ZJ" 


CO  0 
0  CO 


(A-14) 


and  the  modal  Hamiltonian  is 

9i  =  =  -^co(yi+y2)  (A-15) 


Finally,  the  equations  of  variation  are 


39? 

39? 


(A-16) 
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Case  5:  Poincar^  exponents  are  a  pair  of  zeros,  the 
degenerate  case,  and  the  symplectic  eigenvectors  are  real. 
The  initial  and  final  form  of  J  are 


J  = 


0  1 
0  0 


(A-17) 


and  the  resulting  S  is 


S  =  -ZJ  = 


0  0 
0  1 


(A-18) 


and  the  modal  Hamiltonian  91  is 

91  =  ^y^Sy  = 


Finally,  the  equations  of  variation  are 


Yz 


d9l 

_  d9l 


=  Y2 
=  0 


(A-19) 


(A-20) 
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Case  6:  Poincare  exponents  are  a  pair  of  zeros,  the 
degenerate  case,  and  the  symplectic  eigenvectors  are 
imaginary.  The  forms  of  J  are 


J  = 


b  1 
0  0 


J" 


b  0 
1  0 


(A-21) 


and  the  resulting  S' '  is 


S"  =  -ZJ"  = 


-1 

0 


0 

0 


(A-22) 


and  the  modal  Hamiltonian  «  is 

9t  .  ^y's"F  - 


Finally,  the  equations  of  variation  are 


yi 

Yi 


-^-0 


.  m 
■9^ 


=  yi 


{A-23) 


(A-24) 
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Appendix  B:  Fortran  Code 


c 

program  section 
c 

c  ****************  "SURFACE  OF  SECTION"  ********************** 
c  -  creates  surface  of  section  plot  for  the  restricted 
c  three  body  perturbation  problem 

c 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (0-2) 
c 

c  ******************  program  commons  ************************* 

c 

common  /data/  xmu, xmua 

common  /ham/  t, x (20, 4) , f (20, 4) , err (20) , hh, nn,mode 
c 

dimension  rdotv  (4),x(20,4),f(20,4),err(20) 
c 

character*10  filenamel, filename2, filenames 
c 

c  ******************  read  input  data  ************************* 

c 

read(*,*)  nic 
c 

c  nic  =  number  of  initial  conditions 
c 

read(*,*)  xmu 
xmua  =  l.dO  -  xmu 
read(*,*)  hh, tmax 
npts  =  dint (tmax/hh) 
read(*,*)  xjac,  syn 
c 

c  ****************  OPEN  OUTPUT  FILES  ************************* 

c  -  file  1  is  general  output 

c  -  file  2  is  input  for  program  period 

c  -  file  3  is  the  surface  of  section  plotfile 

c 

read(*,*)  filenamel 

open (1,  FILE=f ilenamel, STATUS=' UNKNOWN' ) 
read(*,*)  filename2 

open (2,FILE=filename2, STATUS=' UNKNOWN' ) 
read(*,*)  filenames 

open (3, FILE=filename3, STATUS=' UNKNOWN' ) 
c 

c  *********  repeat  input  values  to  data  FILE  ***************** 

c 

write  (1,*)  'mu=',xmu,'  l-mu=',xmua 
write  (1,*)  'step=',hh,'  maximum  time=',tmax 
write  (1,*)  'number  of  points=' , npts 
write  (1,*)  'Jacobian  constant=' , x jac 
c 

do  600  j  =  l,nic 

)t=0 

read(*,*)  xnot,ynot 
write  (1,*) 

write  (1,*)  'x0=',xnot,'  y0=',ynot 
c 

mode  =  0 
nn  =  4 
nxt  =  0 
t  =  O.dO 
c 

c  GET  ql,pl,q2,p2,AND  HAMILTONIAN;  FOR  GIVEN  x0,y0,  AND  JACOBIAN 

**************** 

ql  =  xnot+xmu 
q2  =  ynot 

xham  =  (xmu*xmua-xjac) /2.d0 
rl  =  dsqrt ( (ql-xmu) * (ql-xmu)  +  q2*q2) 
r2  =  dsqrt ( (ql+xmua) * (ql+xmua)  +  q2*q2) 
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d  =  xham  +  xmua/rl  +  xmu/r2 

g  =  q2/ (ql-xmu) 

a  =  0.5d0* (g*g  +  l.dO) 

b  =  - (g*g*xmu  +  g*q2  +  ql) 

c  =  0 . 5d0*g*g*xmu*xmu  +  g*q2*xmu  -  d 

disc  =  b*b  -  4.d0*a*c 

if(  disc  .It.  O.dO)  then 

write  (1,*)  'stop  10  -  no  roots' 
go  to  600 

endif 

p2  =  (-b+syn*dsqrt (disc) ) / (2 .d0*a) 
pi  =  g*(xmu-p2) 
c 

c  *****************  initial  conditions  ***************** 

c 

X  (1,  1)  =  ql 
x(2,l)  =  pi 
x(3,l)  =  q2 
x(4, 1)  =  p2 
c 

c  *****************  initialize  naming  ******************* 

c 

call  haming(nxt) 
c 

c  ***********  turn  off  second  eom  evaluation  ************ 
c 

nxt  =  -nxt 

if(nxt  .ne.  0)  go  to  499 

write  (1/*)  'stop  99  -  unable  to  start  Haming' 
go  to  600 
499  continue 

c 

c  ******************  integration  loop  ******************* 

c 

write  (1, *) 

write  (1,*)  'npts  value  at  n*period' 
do  500  1  =  l,npts 
c 

Q  *****************  CHECK  FOR  ESCAPE  ******»»♦*♦****♦*♦ 

c 

if (nxt  .eq.  1)  then 

rl  =  dsqrt (  (qlc  -  xmu) * (qlc  -  xmu)  +  q2c*q2c  ) 
if(rl  .gt.  5.d0)  then 

write  (1,*)  'stop  29  -  escaped' 
go  to  600 

endif 

endif 

c 

c  *****************  permute  indices  ******************** 

c 

nm3  =  nm2 
nm2  =  nml 
nml  =  nxt 
c 

c  ********  integrate  orbit,  haming  permutes  nxt  **»»•»** 

c 

call  haming(nxt) 
if(i  .eq.  1)  go  to  500 
c 

c  ***************  CALCULATE  R  DOT  V  ******************** 

c 

qld  =  f (1, nxt) 
q2d  =  f(3,nxt) 

rdotv(nxt)  =  (x  (1,  nxt ) -xmu)  *qld  +  x(3,nxt)*q2d 
c 

c  ********  CHECK  FOR  PERI/APOAPSE  CROSSING  ************* 
c 

if  (rdotv (nxt ) *rdotv (nml) .gt . 0 .dO)  go  to  500 
c 

c  ************  CROSSING  HAS  OCCURRED!!!  ****************** 

c  INTERPOLATE  TO  CROSSING  TIME 

c 

)t=!t+l 

frac  =  -rdotv (nxt) / (  rdotv(nxt)  -  rdotv(nml)  ) 
qlc  =  -frac*x (1, nml)  +  (l.dO  +  frac) *x (1, nxt) 
q2c  =  -frac*x (3, nml)  +  (l.dO  +  frac) ♦x (3, nxt) 
xcross  =  qlc-xmu 

ycross  =  q2c 
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write  (3,5)  xcross, ycross 

write  (1,2)  i, x (1, nxt) , xcross, ycross 

500  continue 

write(3,l)  0.0,  0.0 
600  continue 

write  (2, 3)  xmu,npt3 
write(2,4)  xjac,syn 

1  format  ( lx,  dl6 . 9,  3x,  dl6. 9) 

2  format  (lx,  16,  3  (2x,d24. 17)  ) 

3  format (lx, d24 . 17, 4x, i6) 

4  format (lx, d24 . 17, 4x, d24 . 17) 

5  format  (2x,  2  (2x,  e23 . 16)  ) 

close ( 1 ) 
close (2) 
close (3) 

stop 

end 

include  'rhsl.for' 
include  '  flaming,  for' 
include  'h.for' 
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c  ****  SAMPLE  SECTION  INPUT  ****** 
c  number  of  initial  conditions 
c  value  of  mu 

c  integration  time  step  /  total  integration  time 
c  Jacobian  constant  /  sign  for  conversion  to  Szebehely's  eom 
c  general  output  filename 

c  file  for  information  to  be  passed  to  PERIOD 
c  surface  of  section  plotfile 

c  Jefferys's  x  and  y  initial  conditions (enough  pairs  to  equal  first  number) 

c 

c 

6 

0.00095388d0 

O.OOSdO  500. dO 

3.15d0  -l.dO 

slOout 

plOin 

slOplot 

-0.2d0  O.dO 

-0.3d0  O.dO 

-0.4d0  O.dO 

-O.SdO  O.dO 

-0.6d0  O.dO 

-0.7d0  O.dO 
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c 

program  period 
c 

c  *******.*************♦  "PERIOD"  ************************************ 

c  -  finds  symmetric  periodic  orbits  for  Restricted  problem  by 

c  finding  x  and  period  at  zero  to  force  px  and  y  zero  at  period 

c  and  while  keeping  the  hamiltonian  constant 

c 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 
c 

c  ********************  program  commons  ****************************** 

common  /data/  xmu, xmua 

common  /ham/  t, x (20/ 4) , f (20, 4) , err (20) / hh, nn,mode 
c 

dimension  xreal  (4) ,  ximag  (4)  /cerr(2,l),b(2,2) 
dimension  x(20,4),f(20,4),err(20), xww (2) 
dimension  phi (4,  4) ,  xxx (10) ,  rval (2, 4) , rvec(2, 16) 
c 

complex* 16  val (4) ,  vec (4, 4) ,  ww 
c 

equivalence  (val, rval) 
equivalence  (vec, rvec) 
equivalence  (ww,xww) 
c 

character*10  filenamel, f ilename2 
c 

c  ***************************  read  input  data  ******"******************** 

c 

read(*,*)  xmu,npts 
xmua  =  1 . dO  -  xmu 
read(*,*)  xjac,syn 
read(*,*)  tol,maxit 
read(*,*)  xnot, period 
c 


c 

c 

c 

c 


c 

c 

c 


c 


c 

c 

c 

c 

c 

c 


c 

c 

c 


***************************  OPEN  OUTPUT  FILES  ************************* 

-  file  1  is  general  output 

-  file  2  is  input  to  program  exsyrl 

read(*,*)  filenamel 

open(l,FILE=fil ename 1 , STATUS= ' UNKNOWN' ) 
read(*,*)  filename2 

open (2, FILE=filename2, STATUS=' UNKNOWN'  ) 

*******************  repeat  input  values  to  DATA  FILE  ****************** 


write  (1,*) 
write  (1,*) 
write  (1,*) 
write  (1,*) 


' mu=' , xmu, '  l-mu=' , xmua 
'max  iterations=' ,maxit, ' 
'number  of  points=' ,  npts 
'Jacobian  constant=' , x jac 


tolerance=' , tol 


write  (1,*)  'x0=',xnot,'  period=' , period 

ql  =  xnot+xrau 

write  (1,*)  'ql=',ql 

xham  =  (xmu"xmua-x jac) /2 .dO 

write  (1,*)  'xham=',xham 


************************  begin  iteration  loop 


do  500  iter  =  l,maxlt 

************************  SET  UP  initial  state  ************************* 

pi  =  O.OdO 
q2  =  O.OdO 

rl  =  dabs (  ql  -  xmu  ) 
r2  =  dabs (  ql  +  xmua  ) 

p2  =  ql  +  syn*dsqrt (  ql*ql  +  2.d0*(  xraua/rl 
1  +  xmu/r2  +  xheim  ) ) 

**************************  initial  conditions  *♦**•»»*•**•************* 


x(l,l)  =  ql 
x(2,l)  =  pi 
x(3,l)  =  q2 
x(4,l)  =  p2 
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c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


**************************  CALCULATE  TIMESTEP  ************************** 
hh  =  period/ (dble (npts) ) 


write  {!/*)  ' iteration' , iter 

***********************  INITIALIZE  PHI  MATRIX  ************************* 

do  40  i  =  1,4 


do  41  j  =  1,4 

ij  =  4*i+j 

x(ij,l)  =  i 

41 

continue 

x(5*i,l)  =  l.dO 

40 

continue 

******************  initialize  integration  CONSTANTS  ******************* 


mode  =  1 
nn  =  20 
nxt  =  0 
t  =  0,d0 


c 

c  **************************  initialize  haming  ************************** 

c 

call  haming (nxt) 
c 

if (nxt  .ne.  0)  go  to  1000 

write  (1,*)  'failure  to  initialize  -  stop  99' 
go  to  20 

1000  continue 


c 

c  ***************************  INTEGRATION  LOOP  ************************** 
c 

do  1500  i  =  l,npts 

call  haming (nxt) 

1500  continue 


c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 


*************************  extract  error  vector  ************************ 


cerr(l,l)  =  -x(2,nxt) 

Gerr(2,l)  =  -x(3,nxt) 
write  (1,*)  '  errors' 

write  (1,7)  cerr (1, 1) , cerr (2, 1) 


**********************  CALCULATE  CORRECTION  MATRIX 


1 


dp2dql 

b(l,l) 

b(l,2) 

b(2,l> 

b(2,2) 


(pZ-  xmua* (ql-xmu) / (rl*rl*rl) 

-  xmu* (ql+xmua) / (r2*r2*r2) )  /  <p2-ql) 
x(6,nxt)  +  X (18, nxt) *dp2dql 
f (2, nxt) 

x(7,nxt)  +  X (19, nxt) *dp2dql 
f (3, nxt) 


**********************  CALCULATE 


STATE  CORRECTIONS  ******************** 


call  leqt2f (b, 1, 2, 2, cerr, idig, xxx, ier) 


**************************  add  in  CORRECTIONS 


************************* 


write  (i/*)  '  corrections' 

write  (1,7)  cerr  (1, 1) ,  cerr  (2, 1) 
ql  =  ql  +  cerr (1,1) 
period  =  period  +  cerr (2,1) 
write  (1,*)  '  ql=',ql 

write  (1,*)  '  per iod=' ,  period 


************************  CHECK  FOR  CONVERGENCE 


iend  =  0 

if (dabs (cerr (1, 1) )  .gt.  tol)  iend  =  1 
if (dabs (cerr (2, 1) )  .gt.  tol)  iend  =  1 
if (iend  .eq.  0)  go  to  2000 

**********  maximum  iterations  exceeded  without  convergence  ************ 
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500  continue 

write  (1/*)  'Iteration  Limit  Exceeded  -  atop  15' 
go  to  20 


c  ************************  CONVERGED  PROCESSING  ************************* 
c 

2000  continue 

rl  =  dabs (  ql  -  xmu  ) 
r2  =  dabs (  ql  +  xmua  ) 

p2  =  ql  +  syn*d3qrt (  ql*ql  +  2.d0*(  xmua/rl 
1  +  xmu/r2  +  xham  ) ) 

'PROGRAM  CONVERGED  IN',  iter,'  ITERATIONS' 


write (1, *) 
write  (1,  *) 
write (1, *) 
write  (1,  *) 
write  (1,  *) 
write (1, *) 


ql='  ,ql 
pl=',pl 
q2=',q2 

p2=',p2 

period=' , period 


pldot  =  -1 .do* (-1 . d0*p2  +  xmua* (ql-xmu) / (rl*rl*rl) 
1  +  xmu* (ql+xmua) / (r2*r2*r2) ) 

qldot  =  pi  +  q2 

p2dot  =  -l.d0*(pl  +  xmua*q2/ (rl*rl*rl) 

+  xmu*q2/ (r2*r2*r2) ) 
p2  -  ql 


q2dot 
write (1, *) 
write (1, *) 
write (1, *) 
write (1, *) 


qldot=' ,  qldot 
pldot=' , pldot 
q2dot=' , q2dot 
p2dot=' , p2dot 


c 

c 

writed,*)  'phi' 
do  59  i  =  5,20 

writed,*)  x(i,nxt) 

59  continue 

do  60  i  =  1,4 

do  60  j  =  1,4 

phi(j,i)  =  x(4*i  +  j,nxt) 

60  continue 
writed,*)  'PHI' 
do  61  i  =  1,4 

writed,  4)  phi  (i,  1),  phi  (i,  2),  phi  (i,  3),  phi  (i,  4) 

61  continue 


c 

c 

c 

c 

c 

c 


c 

c 

c 


***************  COMPUTE  EIGEN  VALUES  AND  VECTORS  OF  PHI  *************** 


call  devcrg(4,phi, 4, val, vec, 4) 


********************  COMPUTE  POINCARE  EXPONENTS 


do  80  i  =  1,4 

ww  =  val (i) 

complex  log  of  eigenvalue  over  period 

xreal(i)  =  dlog (dsqrt (xww (1) *xww (1 )  +  xww (2) *xww (2) ) )  /period 
ximag(i)  =  datan2  (  xww (2),  xwwd)  )  /period 
80  continue 

write  (1,*)  'EVALUEs' 
do  100  1  =  1,4 

writed, 5)  i,  rval  (1,  i) ,  rval  (2,  i) 

100  continue 

writed,*)  'POINCARE  EXPONENTS' 

writed,*)  '  xreal(i)','  ximag(i)' 

do  120  i  =  1,4 

writed, 5)  1,  xreal  d) ,  xlmag  (i) 

120  continue 

write  d,*)  'EVECTORS' 
do  140  i  =  1,16 

if(  d-l)/4  .eq.  d+2) /4)  then 

writed, 5)  d  +  3)/4  ,  rvec  (1,  i) ,  rvec(2,  i) 

else 

writed, 6)  rvec  (1,  i) ,  rvec  (2,  i) 

endif 

140  continue 


c 

c  ************  OUTPUT  DATA  FOR  NEXT  PROGRAM  ************************** 
write  (2,  8)  xmu,xjac 
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write (2, 3)  period/ npts 
write(2,2)  ql,pl,q2,p2 
write (2/2)  qldot/ pldot/ q2dot/ p2dot 
do  160  i  =  1/4 

write  (2/  2)  phi  (i  / 1)  /phi  (i  /  2)  /phi  (i/  3)  /phi  (i/  4) 
160  continue 

do  180  i  =  1/4 

write (2/ 6)  rval (1, i) / rval <2/ i) 

180  continue 

do  220  i  =  1/  16 

write(2/6)  rvec  (1/ i )  /  rvec  (2/ i) 

220  continue 

2  format (IX/ 4 (d24 . 17/ 2x)  ) 

3  format  (Ix/ d24 . 17/ 4x/ i6) 

4  format (IX/ 4 (dl8.ll/ 2x) ) 

5  format (IX/ il/ 2X/ 2 (d24 . 17/ 2x) ) 

6  format (4x/ 2 (d24 . 17/ 2x) ) 

7  format (16x/ 2 (d24 . 17/ 2x) ) 

8  format (Ix/ 2 (d24 . 17/ 2x) ) 

20  continue 

close (1) 
close (2) 

stop 

end 

include  'rhsl.for' 
include  'haminq.for' 
include  'h.for' 
include  'leqt2f.for' 
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c  **************  SAMPLE  PERIOD  INPUT  *************** 

c  mu  /  number  of  timesteps (most  likely  reduced  from  SECTION  output) 
c  Jacobian  constant  /  conversion  sign 

c  iteration  tolerance  /  maximum  number  of  iterat ions 

c  initial  guess  for  x  /  initial  guess  for  pert -■d 

c  general  output 

c  file  for  data  to  be  passed  to  program  EXSYRL 

c 

c 


0.95388000000U00005D-03 


4000 


0 . 315D+01  -0 . lOOOOOOOOOOOOOOOOD+01 

l.d-10  20 

-0.51366d0  6.25658d0 


plOout 

exlOin 
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c 

program  exsyrl 
c 

c  *****************  "EXTENDED,  SYMPLECTIC,  REAL"  ****************** 
c  -  program  to  calculate  the  extended  eigenvector 
c  for  a  repeated  root,  symplectitize  the  eigenvector 
c  matrix,  and  create  the  equivalent  real  eigenvector  and 
c  J  matrices  from  the  symplectic 

c 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 
c 

c  ********************  program  commons  ************************** 

C 

dimension  xww (2) , xwc (2) , itype (2) 
dimension  phi (4,4), rval (2, 4) , rvec (2, 16) 
c 

complex*  16  xxx  (10),a(2,2),b(2,l),  ww,  wc,  alpha  (2) ,  error 
complex*  16  val  <4)  ,  fvec  (4, 4)  ,  x  j  (4,  4) ,  evec(4,  4) ,  tzvec  (4,4) 
complex*16  evec2 (4, 4) , x j2 (4, 4) 
c 

equivalence  (val, rval) 
equivalence  (fvec, rvec) 
equivalence  (ww, xww) 
equivalence  (wc,xwc) 
c 

character*10  filenamel 
character*10  filename2 
c 

C  ***************************  read  input  data  *************************** 

c  -  data  should  be  manually  rearranged  so  any  degenerate  mode 

c  eigenvalues  and  eigenvectors  are  last 

c 

read(*,*)  xmu,xjac 

read(*,*)  period,  npts 

read(*,‘)  ql,pl,q2,p2 

read(», *)  qldot, pldot, q2dot, p2dot 

do  20  i  =  1,4 

read  (*,  *)  phi  (1,1),  phi  (i,  2) ,  phi  (i ,  3)  ,phi  (i,  4) 

20  continue 

do  40  i  =  1,4 

read(*,*)  rval  (1,  i) ,  rval  :2,  i) 

40  continue 

do  60  i  =  1,16 

read(*,*)  rvec (1,  i) ,  rvec (2,  i) 

60  continue 

read(*,»)  itype(l),  itype(2) 
c  ********  type  0  =  degenerate 

c  type  1  =  positive/negative  real 

c  type  2  =  positive/negative  imaginary 

c 

c  ***********************  open  output  files  ************************** 

c  -  file  1  is  general  output 

c  -  file  2  is  input  to  programs  floquet  and  hamiltonian 

c 

read(*,*)  filenamel 

open (l,FILE=fi lenamel,STATUS=' UNKNOWN'  ) 
read(*,*)  fllename2 

open (2, FILE=f ilename2, STAT0S=' UNKNOWN'  ) 
c 

Q  ********************  repeat  input  to  output  ♦♦*♦*♦*♦*■*♦*♦****♦*♦*♦♦♦ 

c 

write (1,*)  '  period=' , period 

writed,*)  '  ql=',ql 

c 

c  **************  find  POINCARE  EXPONENTS  AND  J  MATRIX  ************** 
c 

write (1, *) 

writed,*) 'POINCARE  EXPONENTS' 

writed,*)  '  REAL  IMAJ' 

do  80  i  =  1,4 

do  100  j  =  1,4 

xj(j,l)  =  (O.dO,  O.dO) 

100  continue 

ww  =  val  ( 1 ) 
c 
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c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


complex  log  of  eigenvalue  over  period 

iii  =  int ( (i+1) /2) 
if  (itype(iii)  .eq.  0)  then 
xwc (if  =  0 . dO 
xwc(2)  =  O.dO 

elseif  (itype(iii)  .eq.  1)  then 

xwc(l)  =  dlog (dsqrt (xww (1 ) *xww (1) 

S  +  xww (2) *xww (2) ) )  /period 

xwc (2)  =  O.dO 

elseif  (itype(iii)  .eq.  2)  then 
xwc (1)  =  O.dO 

xwc (2)  =  datan2 (  xww (2),  xww(l)  )  /period 

else 
endi  f 

X  j  ( i ,  i )  =  wc 
write(l,2)  i,Mc 
80  continue 

if  (itype(2)  .eq.  0)  x j (3, 4) = (1 . dO,  O.dO) 

*****  assumes  no  more  than  one  degenerate  node  -  pair  placed  second 
write (1, *) 

write (1,  *) ' J  MATRIX' 

write (1,*)  '  REAl  IMAJ' 

do  120  i  =  1,4 

do  120  j  =  1,4 

if{j  eq.  1)  then 

write(l,2)  i,xj(3,i) 

e-  ue 

write(l,3)  xj(j,i) 

endif 

120  continue 

************  FIND  EXTENDED  EIGENVECTOR  IF  THERE  IS  ONE  ************ 

if (i type (2)  .ne.  0)  goto  140 

a(l,l)  =  phi(l,l)  -  val(4) 

a(l,2)  =  phi(l,4) 

a(2,l)  =  phi(3,l) 

a{2,2)  =  phi(3,4) 

b(l,l)  =  fvec  (1,  3)  *period 

b(2,l)  =  fvec (3,  3)  *period 

call  cleqt2f (a, 1, 2, 2, b, idig, xxx, ier) 

if  (ier  .eq.  0)  go  to  160 

write  (1,*)  'stop  129  -  MATRIX  IS  SINGULAR' 
goto  2000 
160  continue 

f''ec(l,4)  =  b(l,l) 
f-.ec(2,4)  =  (O.dO,  O.dO) 
fvec(3,4)  =  (O.dO,  O.dO) 
fvec(4,4)  =  b(2,l) 

140  write (1, *) 

write (1,*)' EIGENVECTOR  COLUMNS  AFTER  EXTENDED  FOUND  (IF  NEEDED)' 
write  (1,*)  '  REAL  IMAJ' 

do  180  i  =  1,4 

do  180  j  =  1,4 

if(j  .eq.  1)  then 

write(l,2)  i,fvec(j,i) 

else 

write(l,3)  fvec(j,i) 

endif 

180  continue 

**********  CHECK  TO  SEE  IF  EVECTORS  AND  EVALUES  ARE  VALID  ******* 

call  chec)t  (phi,  val,  fvec,  period,  error) 
write ( 1, * ) 

write  (1,*)  'max  error  in  evalue  and  evector  chec)c' 
writed,*)  ' error=' , error 

******************  FIND  SYMPLECTIC  EIGENVECTORS  ****************** 

call  symplec (fvec, tzvec, error ) 
write (1, *) 

writed,*)  'max  error  in  symplectic  F''T*Z*F' 
writed,*)  '  error=' ,  error 
ww  =  error 


110 


o  o 


if  (dabs (xww ( 1) )  .gt.  l.d-10)  goto  200 
if  (dabs (xww (2) )  .gt.  l.d-10)  goto  200 
goto  220 
c 

c  **********  get  multiplication  factors  and  correct  f  matrix  ****** 
c 

200  alpha(l)  =  (l.dO,  0 .dO) /cdsqrt (tzvec  (1/ 2) ) 
alpha(2)  =  (l.dO,  0 .  dO) /cdsqrt  (tzvec  (3,  4) ) 
writed,*)  '  alphal=' ,  alpha  (1 ) 
writed,*)  '  alpha2=' ,  alpha  (2) 
do  240  i  =1,4 

evecd/l)  =  fvec  (i,  1)  *alpha  (1) 
evec(l,2)  =  fvec  (i,  2) ‘alpha  (1) 
evec(i,3)  =  fvec  d,  3) ‘alpha  (2) 
evec{i,4)  =  fvec  (i,  4) ‘alpha  (2) 

240  continue 

call  symplec (evec, tzvec, error) 
write (1, ‘) 

write  (1,‘)  'max  error  in  symplectic  E'T‘Z‘E' 
write  (1,‘)  '  error=' ,  error 
goto  260 

220  do  280  i  =  1,4 

do  280  j  =  1,4 

evec(j,i)  =  fvec(j,i) 

280  continue 
260  writed, ‘) 

write d,‘)' symplectic  EIGENVECTOR  COLUMNS' 
write  (1,‘)  '  REAL  IMAJ' 

do  300  i  =  1,4 

do  300  j  =  1,4 

lf(j  .eq.  1)  then 

writed, 2)  i,evec(j,i) 

else 

writed, 3)  evec(j,i) 

endi  f 

300  continue 

*****************  FIND  real  E  and  j  MATRICES  ****************** 

call  real (evec, x j, evec2, x j2, itype) 

c  *******************  CHECK  IF  STILL  SYMPLECTIC  ****************** 
c 

call  symplec (evec2, tzvec, error) 
write (1, ‘) 

write  d,‘)  'max  error  in  symplectic  after  E  made  real' 
writed, ‘)  '  error=' ,  error 
write (1, *) 

write (1,‘)' SYMPLECTIC  REAL  EIGENVECTOR  COLUMNS' 
write  (1,‘)  '  REAL  IMAJ' 

do  320  i  =  1,4 

do  320  j  =  1, 4 

if ( j  . eq.  1 )  then 

writed, 2)  i,evec2(j,t) 

else 

write  (1,3)  evec2(3,i) 

endif 

320  continue 
write (1, ‘) 

write (1,‘)' REAL  J  MATRIX  COLUMNS' 

write (1,‘)  '  REAL  IMAJ' 

do  340  i  =  1,4 

do  340  j  =  1,4 

i f ( j  .eq.  1 )  then 

writed, 2)  i,xj2(j,i) 

else 

writed, 3)  xj2(3,i) 

endif 

340  continue 
c 

c  ******************.*  OUTPUT  FOR  NEXT  PROGRAM  ****************** 

c 

write (2, 3)  xmu,xjac 
write (2, 6)  period, npts 
write(2,7)  ql,pl,q2,p2 
write(2,8)  itype  (1)  ,  itype  (2) 
do  400  1  =  1,4 

do  400  j  =  1,4 
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write(2,'7)  evec2(j,i) 

400  continue 

do  420  i  =  1,4 

do  420  j  =  1,4 

write(2,’7)  xj2(j,l) 

420  continue 
2000  continue 


2  format  ( lx,  il,  2  (2x,  d24  .  m  ) 

3  format (2x, 2 (2x, d24 . 17) ) 

5  format (lx, d24 . 17) 

6  format (lx, d24 . 17, 4x, i6) 

7  format (lx, 4 (2x, d24 . 17) ) 

8  format (lx, 12, 2x, 12) 

close (1) 
close (2) 

stop 

end 

include  ' cleqt2f . for' 
include  ' symplectic. for' 
include  'chec)c.for' 
include  'real. for' 
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C  *********  SAMPLE  EXSYRL  INPUT  ********** ‘ *** 
c  mu  /  Jacobian  constant 
c  period  /  number  of  timesteps 
c  initial  values  of  ql,  pi,  q2,  p2 
c  initial  values  of  qldot,  pldot,  q2dot,  p2dot 
c  phi  by  row 

c  real  &  imaginary  components  of  eigenvalues 

c  real  &  imaginary  components  of  eigenvectors (matched  in  order  with  eigenvalues) 

c  types  of  Poicare  exponent  pairs 

c  general  output  file 

c  file  for  data  to  be  passed  to  FLOQUET  and  HAMILTONIAN 
c 


0. 

0. 

-0, 

0. 

0. 

0. 

-0, 

0. 


95388000000000005D-03 
62565804 1 1351884 6D+01 


0.315D+01 
4000 

O.OOOOOOOOOOOOOD+00 
0.2265881088747D+01 
0.1755076922110D-01 
0.8609790734677D+00 
0.5556742447948D-01 
0.3961686779099D-01 
98640988974257526D+00  -0 . 1 6430316317961827D+00 

98640988974257526D+00  0. 1 643031631 7961827D+00 

OD+00  O.OOOOOOOOOOOOOOOOOD+00 
OD+00  O.OOOOOOOOOOOOOOOOOD+00 
70221866516062548D-12  -0 . 10442224677270372D-01 

0 
0 
0 


5127081520827D+00 
OOOOOOOOOOOOOD+00 
8609790735321D+00 
1686687698163D+03 
6775180018329D+02 
3138081069484D+00 
0 
0 
1 
1 

-0 
0 

-0 
0 

-0 
0 


O.OOOOOOOOOOOOOD+00 

-0.1003813736061D+01 

-0.3961686779363D-01 

-0.3138081071241D+00 

0.1125430816271D+01 

0.8942606411584D-01 


lOOOOOOOOOOOOOOOOD+01 
39970546784512651D+00 
17951039960051318D-11 
70221866516062548D-12 
lOOOOOOOOOOOOOOOOD+01 
-0.39970546784512651D+00 
0. 17951039960051318D-11 
O.OOOOOOOOOOOOOOOOOD+00 
0.22658810887474097D+01 
-0.10038137360616357D+01 
O.OOOOOOOOOOOOOOOOOD+00 
-0.17366902140924529D-05 
O.lOOOOOOOOOOOOOOOOD+01 
-0.44301253975986643D+00 
0.43449248273103301D-05 
2  0 


OOOOOOOOOOOOOOOOOD+00 
76641470947436119D-12 
23570946055484081D-01 
0.104 42224 677270372D-01 
O.OOOOOOOOOOOOOOOOOD+00 
0.76641470947436119D-12 
-0.23570946055484081D-01 
O.dO 
O.dO 
O.dO 
O.dO 

O.OOOOOOOOOOOOOOOOOD+00 

O.OOOOOOOOOOOOOOOOOD+OO 

O.OOOOOOOOOOOOOOOOOD+00 

O.OOOOOOOOOOOOOOOOOD+OO 


exlOout 

flOin 


-0. 1516521888144D+01 
O.OOOOOOOOOOOOOD+00 
-0.5556742445356D-01 
0.6775180018322D+02 
-0.2722871816792D+02 
0.1125430816198D+01 
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c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 


c 

o 

c 

c 


c 

c 

c 


program  floquet 

**********************  "FLOQUET"  ************************** 

-  turns  state  elements  and  elements  of  E  into  sets 
of  100  Fourier  coefficients  in  order  to  create 
exact  solution 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 

********************  program  commons  ********************** 

common  /data/  xmu, xmua 

common  /ham/  t, x (20, 4) , f (20, 4) / err (20) / hh, nn,mode 
common  /fdat/  xj(4,4) 

dimension  x(20,4),f(20,4),err(20),  xww  (2) 
dimension  s(4,200),v(4,4,200),  itype  (2) 
dimension  clc  (101) ,  sk  (101) ,  fxn  (200) 
dimension  vec  (4,4),  vecf  (4,4),diff(4,4),xji4,4) 

character*10  filenamel,  filenaire2 

********************  read  input  data  ******************** 

read(*,*)  xmu,yjac 
xmua  =  1 . dO  -  xmu 
read(*,*)  period,  npts 
hh  =  period/ (dble (npts) ) 
itrip  =  npts/200 
read(*,*)  qlo,  plo,  q2o,  p2o 
read(*,*)  itype  (1) ,  itype  (2) 

******  itest  makes  certain  all  elements  of  E  and  J  are  real 

itest  =  0 
do  40  i  =  1,4 

do  40  j  =  1,4 

read(*,*)  xww (1) , xww (2) 
if  (xww(2)  .gt.  l.d-10)  itest  =  1 
vec  ( j, i)  =  xww (1) 

40  continue 

do  60  i  =  1,4 

do  60  j  =  1,4 

read(*,*)  xww (1) ,  xww (2) 
if  (xww(2)  -gt.  l.d-10)  itest  =  1 
X]  ( j,i)  =  xww(l) 

60  continue 

if  (itest  .eq.  0)  goto  65 

writed,*)  'There  is  an  imaginary  component  of  E  or  J' 
write  (1,*)  'stop' 
goto  2000 
65  continue 


c 

c 

c 

c 

c 


c 

c 

c 


-  file  1  is  general  output 

-  file  2  is  input  to  program  exact 


read(*,*)  filenamel 

open(l,FILE=filenamel, STATUS=' UNKNOWN' ) 
read(*,*)  filename2 

open (2, FILE=filename2, STATUS=' UNKNOWN'  ) 

*******************  repeat  input  VALUES  TO  DATA  FILE  ****************** 


write (1, *) 

'  xmu=' , xmu 

write (1, *) 

'  xjac=',xjac 

write (1, *) 

'  period=' , period 

write (1, ♦) 

'  points=' , npts 

write (1, *) 

'  trip=', itrip 

write  (1, *) 

' timestep=' ,  hh 

write  (1, *) 

'initial  conditions 

write (1, *) 

'  qlo=',qlo 

write  (1, *) 

'  plo=',plo 

write (1, *) 

'  q2o=',q2o 

write (1, *) 

'  p2o=',p2o 
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E  MATRIX  BY  COLUMN' 


C 

c 

c 


c 

c 

c 


c 

c 

o 


c 

c 

c 


c 

c 

c 


write (1, *)  ' 
do  80  i  =  1,4 

do  80  j  =  1,4 

if(j  .eq.  1)  then 

write(l,2)  i,vec(j,i) 

else 

write(l,3)  vec(j,i) 

endif 

80  continue 

writed,*)  '  J  MATRIX  BY  COLUMN' 

do  100  i  =  1,4 

do  100  j  =  1, 4 

if ( j  .eq.  1)  then 

writed, 2)  i,xj(j,i) 

else 

writed, 3)  xj(j,i) 

endif 

100  continue 

*ir-k1r-ii***-k*-IHi*********1r  SET  UP  INITIAL  STATE  *************** 

X (1, 1)  =  qlo 
x(2,l)  =  plo 
X (3, 1 )  =  q2o 
X  (4, 1)  =  p2o 

**************  INITIALIZE  E  AT  t  (o)  ********************** 

do  120  i  =  1,4 

do  120  j  =  1,4 

x(i*4+j,l)  =  vec(j,i) 

120  continue 
mode  =  1 
nn  =  20 
nxt  =  0 
t  =  O.dO 

*****************  initialize  HAMING  ********************** 

call  haming(nxt) 
if (nxt  .ne.  0)  goto  140 

writed,*)  'FAILURE  TO  INITIALIZE  -  STOP  99' 
writed, 7)  fd,l),f(2,l) 
writed, 7)  f  (3,  1) ,  f  (4, 1) 
goto  2000 
140  continue 

******************  begin  INTEGRATION  LOOP  ****************** 

do  160  i  =  1,200 

do  180  j  =  1,4 

s(j,i)  =  x(j,nxt) 
do  180  )^  =  1,4 

v(l^,j,i)  =  x(  (j*4)+)c,nxt) 

180  continue 

do  200  ii  =  l,itrlp 
call  haraing(nxt) 

200  continue 

160  continue 

**************  feed  STATE/EVECTORS  to  FOURIER  ************ 

write  (2,  5)  xmu 

write (2, 6)  period, npts 

write  (2,  8)  itype  (1 ) ,  itype  (2) 

xnot  =  qlo  -  x:nu 

ynot  =  q2o 

write  (2, 7)  xnot, ynot 

write (2,*)  'insert  increase/decrease  to  xnot,xjac' 

write(2,7)  xjac,  -l.dO 

write (2,*)  'insert  itrip  value  here' 

write (2,*)  'insert  two  filenames  here' 

do  220  i  =  1,4 

do  240  j  =  1,200 

fxn ( j)  =  s (i, j) 

240  continue 

call  fourier  (fxn,  c)t,  s)c,  100) 
do  260  j  =  1, 100 
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write(2,7)  ck(j),s)c(j) 

260  continue 

220  continue 

do  280  i  =  1,4 

do  280  j  =  1,4 

do  300  k  =  1,200 

fxn(k)  =  v(  j,i,k) 

300  continue 

call  fourier (fxn, ck, sk, 100) 
do  320  k  =  1,100 

wrlte(2,7)  ck(k),3k(k) 

320  continue 

280  continue 
c 

c  **************  extract  final  state  ******************* 

c 

qlt  =  x(l,nxt) 
pit  =  x(2,nxt) 
q2t  =  x(3,nxt) 
p2t  =  x(4,nxt) 
do  340  i  =  1,4 

do  340  j  =  1,4 

vecf{],i)  =  x(i*4+j,nxt) 

340  continue 


c 

c 

c 


c 

c 

c 


*******************  final  state  conditions 


360 


writed,*)  'STATE 

AT  Tf' 

writed,*)  ' 

qlt=' ,qlt 

writed,*)  ' 

plt=' ,plt 

writed,*)  ' 

q2t=',q2t 

writed,*)  ' 

p2t=',p2t 

writed,*)  ' 

E  MATRIX  BY  COLUMN 

do  360  i  =  1,4 

do  360  j  =  1, 

4 

if(j  .eq.  1)  then 

write(l,2)  i,vecf(j,i) 

else 

write  (1,3)  vecf(3,i) 

endif 

diff(j,i)  =  vecf  ( j,  l)-vec(  j,  i) 

continue 


********  difference  in  initial  and  final  conditions  ************** 


write (1, *) 
write (1, *) 

writed,*)  '  REAL  E(t)-E(0)' 

do  400  i  =  1, 4 

do  400  j  =  1,4 

if(j  .eq.  1)  then 

writed, 2)  i,diff(j,i) 

else 

writed, 3)  diff(j,i) 

endif 

400  continue 


2  format (lx, il,2x,d24. 17) 

3  format (4x, d24 . 17) 

5  format (2x, d24 . 17) 

6  format  (2x,  d24 . 17,  2x,  i6) 

7  format (2x, 2 (2x, d24 . 17) ) 

8  format (2x, i2, 4x, 12) 

2000  continue 

close (1) 
close  (2) 

stop 

end 

Include  'rhs2.for' 
include  'haming.for' 
include  'h.for' 
include  ' fourier . for' 
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c  *********  SAMPLE  FLOQUET  INPUT  ************ 
c  mu  /  Jacobian  constant 
c  period  /  number  of  integration  steps 
c  initial  values  of  ql,  pi,  q2,  p2 
c  types  of  Poincare  exponent  pairs 

c  E  matrix  by  column  (real  &  imaginary  components  -  imaginary  should  be  0.0) 
c  J  matrix  by  column  (real  &  imaginary  components  -  imaginary  should  be  0.0) 
c  general  output  file 
c  file  for  passing  data  to  EXACT 
c 


0.953880000000000050-0. 
0. 625658041135188460+01 
-0.512708152082727580+00 
-0.151652188814436340+01 
2  -1 

0.326832409482314960+00 
0.000000000000000000+00 
0.239879227592609820-10 
-0.737749792901919550+00 
-0.219788076627480680-10 
0.312991167671130730+02 
-0.125104281105381750+02 
0.561851676295077600-10 
0.831333063658177270+01 
O.OOQOOOOOOQOQOOOOQD+00 
0.000000000000000000+00 
-0.207986412632198780+02 
0.000000000000000000+00 
-0.111021275865798170+01 
0.491838173956355040+00 
0.000000000000000000+00 
0.000000000000000000+00 
0.263804700727756570-01 
0.000000000000000000+00 
0.000000000000000000+00 
-0.263804700727756570-01 
0.000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0.100000000000000000+01 
0.000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
f lOouta 
fq251n 


0.3150+01 

4000 

0.000000000000000000+00 


0.277555756156289140-16 
0.000000000000000000+00 
0.000000000000000000+00 
-0.111022302462515650-15 
0.000000000000000000+00 
0.355271367880050090-14 
-0.888178419700125230-15 
0.000000000000000000+00 
0.000000000000000000+00 
0 . 000000000000000000+00 
0.000000000000000000+00 
0. 000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0. 000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0 . 000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0. 000000000000000000+00 
O.OOOOOOOOOOOOOOOOOD+00 
0.000000000000000000+00 
0.000000000000000000+00 
0. 000000000000000000+00 
0. 000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 
0.000000000000000000+00 


0. 000000000000000000+00 
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c 

c 

c 

c 

c 

c 


c 

c 


c 


c 


c 

c 

c 

c 


c 

o 

c 


c 

c 

c 

c 

c 


c 

c 

c 


program  hamiltonian 

********************  "HAMILTONIAN"  ************************ 

-  turns  periodic  coefficients  needed  to  expand  hamiltonian 
into  sets  of  100  fourier  coefficients 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 

*********»*.***  program  commons  ************************ 

common  /data/  xmu,xmua 

common  /ham/  t, x (20, 4) , f (20, 4) , err (20) , hh, nn, mode 
common  /fdat/  xj(4,4) 

dimension  x(20,4),f(20,4),err(20), xww (2) 
dimension  c)t(101),s)c(101),fxn{200) 
dimension  c(20,200),s(4),  itype  (2) ,  w  (2) 
dimension  vec  (4,  4) ,  vecf  (4,4),diff(4,4),xj(4,4) 
dimension  e  (4,  4) ,  xh3  (4,4,4),tc(4,4,4) 

equivalence  (ww, xww) 
character*10  filenamel, f ilename2 


********************  read  input  data  ******************** 

read (*, *)  xmu 

xmua  =  1 . dO  -  xmu 

read(*,*)  period,  npts 

itrip  =  npts/200 

hh  =  period/ (dble (npts) ) 

read(*,*)  qlo, plo,  q2o,p2o 

read(*,*)  itype  (1) ,  itype  (2) 

******  itest  maltes  certain  all  elements  of  E  and  J  are  real 

itest  =  0 
do  40  i  =  1,4 

do  40  j  =  1,4 

read ( * , * )  xww ( 1 ) , xww ( 2 ) 
if  (xww(2)  .gt.  l.d-10)  itest  =  1 
vec(j,i)  =  xww(l) 

40  continue 

do  60  i  =  1,4 

do  60  j  =  1,4 

read(*,*)  xww  (1) ,  xww  (2) 
if  (xww(2)  .gt.  l.d-10)  itest  =  1 
xj(  j,i)  =  xww(l) 

60  continue 

if  (itest  .eq.  0)  goto  65 

writed,*)  'There  is  an  imaginary  component  of  E  or  J' 
write  (1, *)  ' stop' 
goto  2000 
65  continue 

***************************  OPEN  OUTPUT  FILES  ************************* 

-  file  1  is  general  output 

-  file  2  is  input  to  program  expand 

read(*,*)  filenamel 

open (1, FILE=filenamel, STATUS=' UNKNOWN' ) 
read(*,*)  fllename2 

open (2, FILE=filename2, STATUS^' UNKNOWN' ) 

*********"*********  repeat  input  values  TO  DATA  FILE  ****************** 

write (1,*)  '  xmu=' , xmu 

writed,*)  '  xjac=',xjac 

writed,*)  '  per  iod=' ,  period 
writed,*)  '  points=' ,npts 
writed,*)  '  trip=', itrip 
writed,*)  '  timestep=' ,hh 
writed,*)  'initial  conditions' 
writed,*)  '  qlo=',qlo 

writed,*)  '  plo=',plo 

writed,*)  '  q2o=',q2o 
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c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


write (1,*)  '  p2o=',p2o 

write  (1,*)  '  E  MATRIX  BY  COLUMN' 

do  80  i  =  1,4 

do  80  j  =  1,4 

if(j  .eq.  1)  then 

write  (1,2)  i,  vec(  j,  i) 

else 

write(l,3)  vec(j,i) 

endif 

80  continue 

write (1,*)  '  J  MATRIX  BY  COLUMN' 

do  100  i  =  1,4 

do  100  j  =  1, 4 

if(j  .eq.  1)  then 

write  (1,  2)  i,xj  ( j,  i) 

else 

write(l,3)  xj(j,i) 

endif 

100  continue 

**********************  SET  UP  INITIAL  STATE  *************** 

x(l,l)  =  qlo 
X (2, 1)  =  plo 
x(3, 1)  =  q2o 
x(4,l)  =  p2o 

**************  initialize  f  at  t(o)  ********************** 

do  120  i  =  1,4 

do  120  j  =  1,4 

X  <i*4+ j, 1)  =  vec ( j, i) 

120  continue 
mode  =  1 
nn  =  20 
nxt  =  0 
t  =  O.dO 

*****************  initialize  HAMING  ********************** 

call  haming(nxt) 
if (nxt  .ne.  0)  goto  140 

write (1,*)  'FAILURE  TO  INITIALIZE  -  STOP  99' 
writed,-?)  f  (1,1), f  (2,1) 
write(l,7)  f(3,l),f(4,l) 
goto  2000 
140  continue 

******************  begin  integration  loop  ****************** 

do  200  i  =  1,200 

do  220  3  =  1,4 

3(j)  =  x(3,nxt) 
do  220  )t  =  1,4 

e()t,  j)  =  X  (3*4+]t,nxt) 

220  continue 

********************  COMPUTE  THIRD  ORDER  TENSOR  ************** 

do  240  j  =  1,4 

do  240  )t  =  1,4 

do  240  m  =  1,4 

xh3(j,)c,m)  =  h  (s,  3,  j,)c,m,  0,  0) 
tc(j,)c,m)  =  O.dO 

240  continue 

****************  COMPUTE  PERIODIC  COEFFICIENTS  ************** 

-  first  loop  variable  tc  - 

do  260  3  =  1,4 
do  260  )c  =  1,4 
do  260  m  =  1,4 

-  second  loop  variable  xh3  - 

do  260  33  =  1,4 
do  260  )t)c  =  1,4 
do  260  mm  =  1,4 
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& 

260  continue 

c(l,i)  = 
c(2,i)  = 
c(3,i)  = 
c(4,i)  = 
c(5,i)  = 
c(6,i)  = 

& 

c(7,i)  = 

& 

c(8,i)  = 
c(9,l)  = 

& 

c(10,i)  = 
c(ll,i)  = 
c(12,i)  = 
c<13,i)  = 
c(14,i)  = 
c(15,i)  = 

& 

c(16,i)  = 
c(17,i)  = 
c(18,l)  = 
c(19,i)  = 
c(20,i)  = 
do  360  j 
call 

360  continue 

200  continue 


tc(j,)c,m)  =  tc(j,k,in)  + 

xh3  ( j  j,  kk,inm)  *e  ( j  j,  j)  *e  (kk,  k)  *e  (mm,  m) 

tc(l,l,l)/6.d0 

(tc(l,l,2)+tc(l,2,l)+tc(2,  l,l))/6.d0 
(tc(l,l,3)+tc{l,3,l)+tc(3, 1,  1)  )/6.d0 
(tc(l,  1,  4)  +tc(l,4,  l)+tc(4, 1,1) )  /6.d0 
(tc(l,2,2)+tc(2,  l,2)+tc(2,2, 1) )  /6.d0 
(tc(l,2,  3)+tc(l,3,2)+tc(2,l,3)+tc(2,3, 1) 

+tc(3,  l,2)+tc(3,2,l) )  /6.d0 
(tc{l,2,4)+tc(l,4,2)+tc(2,l,4)+tc(2,4,  1) 

+tc(4,  l,2)+tc(4,2, 1)  )/6.d0 
{tc(l,3,3)+tc(3,l,3)+tc(3,3, 1) )  /6.d0 
{tc(l,3,  4)+tc(l,4,3)+tc(3,l,4)+tc(3,4,l) 
+tc(4,l,3)+tc(4,3,l))/6.d0 
(tc(l,  4,  4)  +tc(4, 1,  4)+tc(4,  4, 1) )  /6.d0 
tc(2,2,2)/6.d0 

(tc(2,2,3)+tc(2,3,2)+tc(3,2,2))/6.d0 
{tc{2,2,  4) +tc  (2,4,2) +tc  (4,2,2)  )/6.d0 
(tc(2,3,3)+tc(3,2,3)+tc(3,3,2))/6.d0 
(tc(2,3,  4)+tc(2,4,3)+tc(3,2,4)+tc(3,  4,2) 
+tc(4,2,3)+tc(4,3,2) )  /6.d0 
(tc(2,  4,  4)+tc(4,2,4)+tc(4,4,2) )  /6.d0 
tc(3,3,3)/6.d0 

(tc(3,3,4)+tc(3,  4,3)+tc(4,3,3)  )/6.d0 
(tc{3,4,4)+tc(4,3,4)+tc(4,4,3))/6.d0 
tc  (4,  4,  4)  /6.d0 
:  l,itrlp 
haming (nxt) 


c 

C  *********  COMPUTE  FOURIER  COEFFICIENTS  FROM  PERIODIC  ONES  ****** 


write (2, 6)  period, npts 

write (2,*)  'insert  modal  initial  displacements  here' 
write (2,*)  'insert  itrip  value  here' 
do  380  i  =  1,3,2 

iii  =  int ( (i+1) /2) 
if  (itype(iii)  .eq.  2)  then 
w<iii)  =  xj(i,i+l) 

w(iii)  =  xj(i,i) 

endif 

380  continue 

write(2,7)  w(l),w(2) 

write (2, 8)  itype (1) , itype (2) 

write  (2,*)  'insert  two  filenames  here' 

do  400  i  =  1,20 

do  420  j  =  1,200 

fxn(  j)  =  c(i,  j) 

420  continue 

call  fourier ( fxn, ck, sk, 100) 
do  440  j  =  1,100 

write(2,7)  ck(j),sk(j) 

440  continue 

400  continue 


c 

c 

c 


EXTRACT  FINAL  STATE 


qlt  =  x(l,nxt) 
pit  =  x(2,nxt) 
q2t  =  x(3,nxt) 
p2t  =  x(4,nxt) 
do  480  i  =  1,4 

do  480  j  =  1,4 

vecf(j,i)  =  x(i*4+j,nxt) 

480  continue 


c 

c 

c 


FINAL  STATE  CONDITIONS 


write (1,*)  'STATE 
writed,*)  ' 
writed,*)  ' 
writed,*)  ' 
writed,*)  ' 
writed,*)  ' 
do  500  i  =  1,4 


AT  Tf' 
qlt=',qlt 
plt=' ,plt 
q2t=' ,q2t 
p2t=' , p2t 

E  MATRIX  BY  COLUMN' 
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do  500  j  =  1,4 

if(j  .eq.  1)  then 

write(l,2)  i,vecf(j,i) 

else 

write(l,3)  vecf(j,l) 

endif 

diff(j,i)  =  vec{j,i)-vecf(j,i) 

500  continue 

writed,*)  '  REAL  E(t)-E(0)' 

do  520  i  =  1,4 

do  520  j  =  1,4 

if(j  .eq.  1)  then 

writed, 2)  i,diff<j,i) 

else 

writed, 3)  diff(j,i) 

endif 

520  continue 

2  format (lx, 11, 2x,d24. 17) 

3  format (4x, d24 . 17) 

6  format  (2x,  d24 . 17,  2x,  16) 

7  format (2x, 2 (2x,d24. 17) ) 

8  format (2x, 12, 4x, i2) 

2000  continue 

close  (1) 
close (2) 

stop 

end 

include  'rhs2.for' 
include  'haming.for' 
include  'h.for' 
include  ' fourier . for' 


C  **************  SAMPLE  HAMILTONIAN  INPUT  ***************** 

c  mu  /  Jacobian  constant 
c  period  /  number  of  integration  steps 
c  initial  values  of  ql,  pi,  q2,  p2 
c  types  of  Poincare  exponent  pairs 

c  E  matrix  by  column  (real  &  imaginary  components  -  imaginary  should  be  0.0) 
c  J  matrix  by  column  (real  &  imaginary  components  -  imaginary  should  be  0.0) 
c  general  output  file 
c  file  for  passing  data  to  EXPANDED 
c 
c 

0.95388000000000005D-03  0.315D+01 

0. 62565804113518846D+01  4000 

-0. 5127081520827D+00  0 . OOOOOOOOOOOOOD+00  0 . OOOOOOOOOOOOOD+00  -0 . 151 65218881 44D+01 

2  -1 

0.32683240948231496D+00  0 . 27 755575615628914D-1 6 

0. OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
0.23987922759260982D-10  0 . OOOOOOOOOOOOOOOOOD+00 

-0.73774979290191955D+00  -0. 1 1102230246251 565D-15 

-0. 21 978807 662748068D- 10  0 . OOOOOOOOOOOOOOOOOD+00 

0.31299116767113073D+02  0 . 35527136788005009D-14 

-0. 125104281 10538 175D+02  -0 . 88817841970012523D-15 

0. 561851 67 629507760D- 10  0  .  OOOOOOOOOOOOOOOOOD+00 

0.83133306365817727D+01  0 . OOOOOOOOOOOOOOOOOD+00 

O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
-0. 2079864 126321 9878D+02  0 . OOOOOOOOOOOOOOOOOD+00 

O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
-0.11102127586579817D+01  0 . OOOOOOOOOOOOOOOOOD+00 

0. 49183817395635504D+00  0 . OOOOOOOOOOOOOOOOOD+00 

O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
0.26380470072775657D-01  0 . OOOOOOOOOOOOOOOOOD+00 

O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
-0.26380470072775657D-01  0 . OOOOOOOOOOOOOOOOOD+00 

0. OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
0. OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
0. lOOOOOOOOOOOOOOOOD+Ol  O.OOOOOOOOOOOOOOOOOD+00 
O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
O.OOOOOOOOOOOOOOOOOD+00  O.OOOOOOOOOOOOOOOOOD+00 
f lOoutb 
hm25in 
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c 

program  exact 
c 

c  ****♦***********.***.*»  "EXACT"  ************************* 

c  -  integrates  a  nearly  periodic  orbit,  subtracts 

c  the  periodic  reference,  transforms  the 

c  result  into  modal  variables,  and  creates  a  plot  file 

c  of  the  modal  variables  bl,  b2,  b3,  and  b4 

c 

implicit  double  precision  (a-h) 
implicit  integer  (1-n) 
implicit  double  precision  (o-z) 
c 

c  ********************  program  commons  ********************** 

c 

common  /data/  xmu, xmua 

common  /ham/  t, x (20, 4) , f (20, 4) , err (20) , hh, nn,mode 
c 

dimension  x(20,4),f(20,4),err(20) 
dimension  cf(20),sinn(100), coss (100) 
dimension  c)t(20,100),s)t(20,100),dx(4,l) 
dimension  e (4, 4) , xxx (SO) , Itype (2) 
c 

equivalence  (ww, xww) 
c 

character*10  filenamel, filename2 
c 

c  **************  read  input  data  ******************** 

c 

read(*,*)  xmu 
xmua  =  1 . dO  -  xmu 
read(*,*)  period, npts 
hh  =  period/ (dble (npts) ) 
read(*,*)  itype  (1) ,  itype  (2) 
read(*,*)  xnot,ynot 
read(*,*)  xd,  xjacd 
xnot  =  xnot+xd 
read(*,*)  xjac,  syn 
xjac  =  xjac  +  xjacd 
read(*,*)  itrip 
c 

c  **********************  OPEN  OUTPUT  FILES  ******************** 
c  -  file  1  is  general  output 

c  -  file  2  is  a  plotfile  of  the  exact  solution 

c 

read(*,*)  filenamel 

open (l,FILE=filenamel,STATUS=' UNKNOWN'  ) 
read(*,*)  filename2 

open(2,FILE=fi 1 ename 2 , STATUS=' UNKNOWN'  ) 
c 

do  60  i  =  1,20 

do  60  j  =  1, 100 

read(*,*)  c)c  (i,  j) , s)t  (i,  j) 

60  continue 
c 

C  *************  GET  ql,pl,q2,p2  FOR  GIVEN  x0,y0,  AND  JACOBIAN  *********** 
c 

ql  =  xnot+xmu 
q2  =  ynot 

xham  =  (xmu*xmua-xjac) /2.d0 

rl  =  dsqrt ( (ql-xmu) * (ql-xmu)  +  q2*q2) 

r2  =  dsqrt ( (ql+xmua) * (ql+xmua)  +  q2*q2) 

d  =  xham  +  xmua/rl  +  xmu/r2 

g  =  q2/ (ql-xmu) 

a  =  0.5d0*(g*g  +  l.dO) 

b  =  - (g*g*xmu  +  g*q2  +  ql) 

c  =  0 . 5d0*g*g*xmu*xmu  +  g*q2*xmu  -  d 

disc  =  b*b  -  4.d0*a*c 

lf(  disc  .It.  O.dO)  then 

write  (1,*)  'stop  10  -  no  roots' 
go  to  2000 

endlf 

p2  =  (-b+syn*dsqrt (disc) ) / (2.d0*a) 
pi  =  g*(xmu-p2) 
c 

c  *******************  repeat  input  VALUES  TO  DATA  FILE  ****************** 

c 

write  (1,*)  'mu=',xmu, '  l-mu=',xmua 
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u  u 


c 


c 


c 

c 

c 

c 


wr  J.  te 

(1,*) 

'step=',hh, '  period^' , period 

write 

(1,*) 

'number  of  points=' , npts 

write 

(1,*) 

write 

(1,  *) 

'initial  conditions  (Jefferys)' 

write 

(1,  *) 

'x0=',xnot,'  y0-',ynot 

write 

(1,  *) 

'Jacobian  constant=' , x jac 

write 

(1,  *) 

write 

(1,  *) 

'initial  conditions  (Szebehely)' 

write 

(1,*) 

'ql^',ql,'  pl=',pl 

write 

(1,*) 

'q2=',q2,'  p2=',p2 

write 

(1,  *) 

'Hamiltonian  constant=' , xham 

mode  ■ 

=  0 

nn  = 

4 

nxt  = 

0 

t  =  0 

■  do 

pi  =  dacos ( 

-1 .dO) 

wO  = 

(2.d0* 

pi )  /'period 

write 

(1, 

'  w0=',w0 

**t*t*i,***t*i,t  initial  conditions  ******************** 


X  ( 1, 

1) 

=  qi 

X  (2, 

1) 

=  pi 

X  (3, 

1) 

=  q2 

X  (4, 

1) 

=  P2 

*♦**.*****.♦****»**..».♦**.  initialize  hawing  ************************* 

call  hamlnq(nxt) 

♦****»»*»»»********»  turn  off  second  eom  evaluation  ******************* 


nxt  =  -nxt 

if(nxt  .ne.  0)  go  to  499 

write  (1,*)  'stop  99  -  unable  to  start  Naming' 
go  to  2000 
499  continue 

»*»!,»»**»*****»», t*  integration  loop  ******»*»****»***‘****** 

do  500  i  =  0,250 

*♦****«»*♦♦*  compute  sin (n*theta) ,  cos (n'theta) ,  n=l,50  ♦***♦***»**♦♦< 

cossd)  =  aoos(w0*t) 

sinn(l)  =  dsin(w0*t) 

coss(2)  =  2 .d0*coss ( 1) *coss (I ) -1  .dO 

sinn(2)  =  2 .dO'sinn ( 1 ) *coss <1 ) 

do  520  J  =  3, 100 

coss(j)  =  2 . d0*coS3 ( j-1 ) *coss ( 1 ) -cos? ( j-2) 
sinn(j)  =  2 . d0*sinn ( j- 1 ) *coss ( 1 ) -Sinn ( j-2) 

520  continue 

********  reassemble  periodic  traj  and  eigenvector  matrix  ********** 


540 

560 


do 

540  j  ^  1,20 

cf  (  j)  =  ck  ( 

Jf  1) 

do  540  k  = 

1, 99 

cf  (j) 

=  cf(j) 

+  ck < j, k+1 ) *coss (k) 

+  sk ( j, k+1) *sinn (k) 

continue 

do 

560  i  =  1,4 

dx  (  j,  1 )  =  X 

( j,  labs 

(nxt)  )  -cf  (  j) 

continue 

************  place  eigenvectors  in  4X4  matrix  for  inversion  **•*••** 

do  580  j  =  1,4 

^o  580  k  =  1, 4 

e(k,  j)  =  cf  (j*4  +  k) 

580  continue 

-1 

**♦  CALCULATE  DELTA  b  WITH  DELTA  b  =  EVECTORS  *  DELTA  X  *** 
idig  =  0 

call  Ieqt2t (e, 1, 4, 4,dx, idig, xxx, ler) 
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if  (i  .ne.  0)  goto  620 
write (1, *) 

write (1,*)  'Initial  Modal  Conditions' 

writed,*)  '  bl  -  b4' 

write(l,5)  dx(l,l) 

writed, 5)  dx(2,l) 

write  (1,5)  dx  (3,  1) 

write  (1,5)  dx  (4,  1) 

620  continue 

write  (2,6)  t,dx(l,l),dx(2,l),dx(3,l),dx(4,l) 
do  640  j  =  l,itrip 

call  haming(nxt) 

640  continue 

500  continue 
c 

c  ****************  extract  final  conditions  ***************** 

c 

write (1, *) 
write (1, *) 

writed,*)  'final  state  minus  initial  state  (Szebehely) 

writed,*)  '  ql,pl,q2,p2' 

writed,*)  x(l,nxt)-ql 

writed,*)  x(2,nxt)-pl 

writed,*)  x(3,nxt)-q2 

writed,*)  x(4,nxt)-p2 

2000  continue 

2  format (2x, 2 (2x, d24 . 17) ) 

4  format  (4x,dl2.5,  2x,d24. 17) 

5  format (4x, d24 . 17) 

6  format (2x, 5 (2x, e23 . 16) ) 


close (1) 
close (2) 
close  (3) 
close  (4) 


stop 

end 

include  'rhsl.for' 
include  '  iiaming.  for' 
include  'h.for' 
include  'leqt2f.for' 
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c  **************  SAMPLE  EXACT  INPUT  *************** 
c  mu 

c  period  /  number  of  integration  steps (for  calculation  of  timestep) 
c  types  of  Poincare  exponent  pairs 

c  initial  value  of  x  and  y  (Jefferys  eom  initial  condition?) 
c  amount  of  change  to  x  /  amount  of  change  to  Jacobian  constant 
c  initial  Jacobian  constant  /  conversion  sign 

c  integration  steps  between  sampled  data(250  total  data  samples) 
c  general  output 

c  plot  file  of  time  and  modal  variables 

c  20  sets  of  100  sine  and  cosine  Fourier  representation  pairs 

c 

c 

0. 9538800000000000D-03 
0.62565804113518473D+01  4000 

2  -1 

-0. 513 6620320827474 lD+00  0 . OOOOOOOOOOOOOOOOOD+00 

O.dO  O.dO 
3.15d0  -l.dO 
630 

fqlout 

exaplotl 

0.14272865159998943D-02  0 . OOOOOOOOOOOOOOOOOD+00 

-0.44658422312811241D+00  0 . 50812964 946800321D-13 

0.95570615732216700D-04  0 . 15916670759175133D-12 

-0. 58321 183772268116D-01  -0 . 24715715585266195D-13 

etc. ...  (a  total  of  2000  pairs  of  numbers) 
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c 


program  expanded 


c 

C  ******************»...  "EXPANDED"  **•*•***♦••*“******•* 

c  -  Using  the  periodic  coefficients  made  by  the  program 
c  hamiltonian,  the  eom  for  the  truncated 

c  hamiitoniaii  case  are  integrated.  A  plotfile  of  modal 

c  variables  bl,  b2,b3,  and  b4  is  created, 
c 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 
c 

c  ****************  program  commons  ********************** 

c 

common  /datarhs/  wO, wl, w2, ck (20, 100) , sk (100,  50) ,  itype (2) 
common  /ham/  t,  x  (20,  4) ,  f  (20,  4)  ,err  (20)  ,hh,nn,mode 
c 

dimension  x(20,4),f(20,4),err(20),  itype  (2) 
dimension  ck(20,  100),sk(20,  100)  ,b0(4),b(4) 
c 

character*10  filenamel, filename2 
c 
c 

c  **************  read  input  data  ******************** 

c 

read(*,*)  period,  npts 
hh  =  period/ (dble (npts) ) 
do  40  i  =  1,4 

read(*,*)  b0(i) 

40  continue 

read(*,*)  itrip 
read(*,*)  wl,w2 
read(*,*)  itype  ( 1 ),  itype  (2) 
c 

Q  ******************  OPEN  OUTPUT  FILES  ********************** 

c  -  file  1  is  general  output 

c  -  file  2  is  a  plotfile  for  the  expanded  case 

c 

read(*,*)  filenamel 

open(l,FILE=fi I ename 1 , STATUS=' UNKNOWN' ) 
readi*,*)  fllename2 

open (2, FILE=f ilename2, STATUS=' UNKNOWN' ) 
c 

do  60  i  =  1,20 

do  60  j  =  1, 100 

read(*,*)  ck (i,  j) ,  sk  (i,  J) 

60  continue 


c 


write (1, *) 

'  orbit  period=' , period 

write (1, *) 

'  #  of  points=' , npts 

write (1, *) 

'  timestep=' ,hh 

write (1, *) 

'initial  state  (MODAL)' 

write (1, *) 

'  bl=',b0(l) 

write (1, *) 

'  b2=',b0(2) 

write ( 1, * ) 

'  b3=',b0(3) 

write (1, *) 

'  b4=',b0(4) 

c 

SET  UP  INITIAL  STATE  *** 

c 

do  80  i  = 

1,4 

x(i,l)  =  b0(i) 
80  continue 

mode  =  0 
nn  =  4 
nxt  =  0 
t  =  O.dO 


c 

pi  =  dacoa(-l.dO) 
wO  =  2.d0*pi/period 
c 

c  ****************  initializing  naming  **************** 

c 


call  haming(nxt) 
lf(nxt  .ne.  0)  go  to  199 
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write  (1,*)  'stop  99  -  unable  to  start  Hamlng' 
go  to  2000 
199  continue 


c  ****************  begin  integration  loop  **************** 

c 

do  220  i  =  1,250 

do  240  j  =  1,4 

b(j)  =  x(j,nxt) 

240  continue 

write(2,4)  t,  b  (I) ,  b  (2) ,  b  (3) ,  b  (4) 
do  300  j  =  I,itrip 

call  haming(nxt) 

300  continue 

220  continue 
c 

c  ******************  extract  final  state  ***************** 
c 

do  340  j  =  1,4 

b(j)  =  x(j,nxt) 

340  continue 
write (1, *) 


write (1, *) 
write (1, *) 
write (1, *) 
write (1, *) 
write (1, *) 


final  state  (Modal)' 
bl=',b(l) 
b2=',b(2) 
b3=',b(3) 
b4=',b(4) 


2000  continue 

2  format  (2x,  2  (2x,  d24  .  n)  ) 

3  format (4x, dl2 . 5, 2x, d24 . 17) 

4  format (2x, 5 (2x, e23 . 16) ) 

close (1) 
close (2) 


stop 

end 

include  '  jiaming.  for' 
include  ' rhs3 . for' 
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C  ****************  SAMPLE  EXPANDED  INPUT  ********************* 
c  period  /  number  of  integration  steps (used  to  determine  the  timestep) 
c  four  initial  modal  displacements (found  in  the  general  output  of  EXACT) 
c  integration  steps  between  sampled  data (250  total  data  samples) 
c  value  of  omega  for  Poincare  exponent  pairs  (0.0  in  degenerate  case) 
c  types  of  Poincare  exponent  pairs 
c  general  output 

c  plot  file  of  time  and  modal  variables 

c  20  sets  of  100  sine  and  cosine  Fourier  representation  pairs 
c 


0.62565804113518473D+01  4000 

0.39599531533278684D-12 
-0.19873186154673566D-15 
-0.14793686800523243D-13 
-0. 101075395487 98371D-13 
630 

-0.2638047004D-01  0 . OOOOOOOOOOOD+00 

2  -1 

-0. 2638047004 9156036D-01 


2 

hmlout 


expplotl 

0.54375306107848423D-01 
0. 1138620842200 1964D-02 
0.73473610830684946D-01 
-0.54706262117310002D-04 
-0.93120150715344754D-01 
etc. ...  (a  total  of  2000 


O.OOOOOOOOOOOOOOOOOD+00 
0.18637175136504425D-13 
0.23774386834207739D-10 
-0.13772226414854315D-12 
0.17098303259355240D-10 
pairs  of  numbers) 
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c 

subroutine  haming(nxt) 
c 

c  hcuning  is  an  ordinary  differential  equations  integrator 
c  it  is  a  fourth  order  predictor-corrector  algorithm 

c  which  means  that  it  carries  along  the  last  four 

c  values  of  the  state  vector,  and  extrapolates  these 

c  values  to  obtain  the  next  value  (the  prediction  part) 

c  and  then  corrects  the  extrapolated  value  to  find  a 

c  new  value  for  the  state  vector, 

c 

c  the  value  nxt  in  the  call  specifies  which  of  the  4  values 

c  of  the  state  vector  is  the  "next"  one. 

c  nxt  is  updated  by  haming  automatically,  and  is  zero  on 

c  the  first  call 

c 

c  the  user  supplies  an  external  routine  rhs(nxt)  which 

c  evaluates  the  equations  of  motion 

c 

common  /ham/  x,  y  (20,  4) ,  f  (20,  4) ,  errest  (20)  ,h,  n,  mode 
C  double  precision  x, y, f, errest, h, hh, xo, tol 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 
c 

dimension  y (20, 4) , f (20, 4) , errest (20) 
c 

c  all  of  the  good  stuff  is  in  this  common  block, 

c  X  is  the  independent  variable  (  time  ) 

c  y(6, 4)  is  the  state  vector-  4  copies  of  it,  with  nxt 

c  pointing  at  the  next  one 

c  f(6, 4)  are  the  equations  of  motion,  again  four  copies 

c  a  call  to  rhs(nxt)  updates  an  entry  in  f 

c  errest  is  an  estimate  of  the  truncation  error  -  normally  not 

c  used 

c  n  is  the  number  of  equations  being  integrated  -  6  or  42  here 

c  h  is  the  time  step 

c  mode  is  0  for  just  EOM,  1  for  both  EOM  and  EOV 

c 

tol  =  O.OOOOOOOOOld+00 

c  switch  on  starting  algoritlitci  or  normal  propagation 

if(nxt)  190,10,200 
c 

c  this  is  hamings  starting  algoritlim. . .  .a  predictor  -  corrector 

c  needs  4  values  of  the  state  vector,  and  you  only  have  one-  the 

c  initial  conditions. 

c  haming  uses  a  Picard  iteration  (slow  and  painful)  to  get  the 

c  other  three. 

c  if  it  fails,  nxt  will  still  be  zero  upon  exit,  otherwise 

c  nxt  will  be  1,  and  you  are  all  set  to  go 

c 

10  xo  =  X 

hh  =  h/2.0d+00 
call  rhs (1) 
do  40  1  =  2,4 

X  =  X  +  hh 

do  20  i  =  l,n 

20  y(i,l)  =  y(i,l-l)  +  hh*f(i,l-l) 
call  rhs(l) 

X  =  X  +  hh 
do  30  i  =  l,n 

30  y(l,l)  =  y(l,l-l)  +  h*f(i,l) 

40  call  rhs(l) 
jsw  =  -10 
50  isw  =  1 

do  120  i  =  l,n 

hh  =  y(l,l)  +  h*(  9.0d+00*f  (i,l)  +  19,0d+00*f  (i,  2) 

1  -  5.0d+00*f (i,3)  +  f(i,4)  )  /  24,0d+00 

if(  dabs (  hh  -  y(i,2))  .It.  tol  )  go  to  70 
isw  =  0 

70  y(i,2)  =  hh 

hh  =  y(i,l)  +  h*(  f(i,l)  +  4 . 0d+00*f  (i,  2)  +  f  (i,  3) ) /3 .  Od+00 
if(  dabs  (  hh-y(i,3))  .It.  tol  )  go  to  90 
isw  =  0 

90  y(i,3)  =  hh 

hh  =  y(i,l)  +  h*(  3.0d+00*f  (i,  1)  +  9. 0d+00*f  (i,  2)  +  9. 0d+00*f  (i,  3) 
1  +  3.0d+00*f (1,4)  )  /  8.0d+00 

if(  dabs (hh-y (i, 4) )  .It.  tol  )  go  to  110 
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isw  =  0 

110  y(i,4)  =  hh 
120  continue 
X  =  xo 

do  130  1  =  2,4 

X  =  X  +  h 

130  call  rhs(l) 

if(isw)  140,  140,  150 
140  jsw  =  jsw  +  1 

if(j3W)  50,280,  280 
150  X  =  xo 
isw  =  1 
jsw  =  1 

do  160  i  =  l,n 
160  errest(i)  =  O.OdO 
nxt  =  1 
go  to  280 
190  jsw  =  2 

nxt  =  iabs(nxt) 
c 

c  this  is  hamings  normal  propagation  loop  - 

c 

200  X  =  X  +  h 

npl  =  mod (nxt, 4)  +  1 
go  to  (210, 230) , isw 
c  permute  the  index  nxt  modulo  4 
210  go  to  (270, 270, 210, 220), nxt 
220  isw  =  2 

230  nm2  =  mod (npl, 4)  +  1 
nml  =  mod(nm2,4)  +  1 
npo  =  mod (nml, 4)  +  1 
c 

c  this  is  the  predictor  part 

c 

do  240  i  =  l,n 

f(i,nm2)  =  y(i,npl)  +  4.0d+00*h*(  2. 0d+00*f  (i,  npo)  -  f(i,nml) 

1  +  2.0d+00»f (i,nm2)  )  /  3.0d+00 

240  y(i,npl)  =  f(i,nm2)  -  0 . 925619835d0*errest  (i) 
c 

c  now  the  corrector  -  fix  up  the  extrapolated  state 
c  based  on  the  better  value  of  the  equations  of  motion 
c 

call  rhs(npl) 
do  250  i  =  l,n 

y(i,npl)  =  (  9 . 0d+00*y  (i,  npo)  -  y(i,nm2)  +  3.0d+00*h*(  f(i,npl) 
1  +  2 . 0d+00*f (i, npo)  -  f(i,nml)  )  )  /  8.0d+00 

errestd)  =  f(i,nm2)  -  y(i,npl) 

250  y(i,npl)  =  y(i,npl)  +  0 . 0743801653d0  •  errest(i) 
go  to  (260, 270) , jsw 
260  call  rhs(npl) 

270  nxt  =  npl 
280  return 
end 
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c 

function  h  (x,  iord,  i,  j,  k,  l,m) 
c 

c  restricted  problem  in  canonical  coordinates 

c 

c  state  vector  x  =  (  ql,  pi,  q2,  p2  ) 

c 

dimension  x(4) 
common  /data/  xmu, xmua 
implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 
c 

c  preliminaries 

c 

qa  =  X  ( 1 )  -  xmu 
qb  =  X ( 1 )  +  xmua 
rl  =  (qa*qa  +  x (3) *x (3) ) ** . 5d0 
r2  =  (qb*qb  +  x ( 3) *x (3) ) ** . 5d0 
c 

c  branch  on  order 

c 

jord  =  iord  +  1 

go  to  (1,  1000,  2000,  3000), jord 
c 

Q  ************************************ 

c  **  ** 

c  Order  Zero  ** 

c  **  ** 

c 

1  continue 

h  =  0.5d0*{x(2)*x(2)  +  x(4)*x{4))  +  x(3)*x(2) 
1  -  xmua/rl  -  xinu/r2 

return 

1000  continue 


c 

c 

c 

c 

c 

c 

c 


c 


•**  ** 

**  Order  One  ** 

**  ** 

************************************ 

rl3  =  rl**3.d0 
r23  =  r2**3.d0 

go  to  (1001,  1002,  1003,  1004),  i 


1001  h  =  “X(4)  +  xmua*qa/rl3  +  xmu*qb/r23 
return 

1002  h  =  x(2)  +  x(3) 
return 

1003  h  =  x(2)  +  xmua*x  (3) /rl3  +  xinu*x  (3) /r23 
return 

1004  h  =  x(4)  -  x(l) 
return 


c 

2000  continue 


c 

c 

c 

c 

c 

c 


**  ** 

**  Order  Two  ** 

*  *  *  * 


rl3  =  rl**3.d0 
r23  =  r2»*3.d0 
rl5  =  rl**5.d0 
r25  =  r2**5.d0 


go 

to 

(2001, 

2002, 

2003, 

2004)  ,i 

2001 

go 

to 

(2011, 

2012, 

2013, 

2014), J 

2002 

go 

to 

(2021, 

2022, 

2023, 

2024), J 

2003 

go 

to 

(2031, 

2032, 

2033, 

2034),! 

2004 

go 

to 

(2041, 

2042, 

2043, 

2044), J 

2011  h  =  xmua/rl3  +  xmu/r23  -3.d0*xmua*qa*qa/rl5 
1  -3.d0*xmu*qb*qb/r25 

return 


x(l)*x(4) 
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2012  h  =  O.dO 
return 

2013  h  =  -3.d0*xmua*qa*x(3)/rl5  -  3 .dO*xmu*qb*x <3) /r25 
return 

2014  h  =  -l.dO 
return 

2021  h  =  O.dO 
return 

2022  h  =  l.dO 
return 

2023  h  =  l.dO 
return 

2024  h  =  O.dO 
return 

2031  go  to  2013 

2032  h  =  l.dO 
return 

2033  h  =  -3.d0*xmua*x(3) *x(3) /rl5  -  3 .d0*xmu»x (3) *x (3) /r25 

1  +  xmua/rl3  +  xmu/r23 

return 

2034  h  =  O.dO 
return 

2041  h  =  -l.dO 
return 

2042  h  =  O.dO 
return 

2043  h  =  O.dO 
return 

2044  h  =  l.dO 
return 

3000  continue 


c 

c 

c 

c 

c 

c 

o 


*  *  ** 

**  Order  Three 

**  *♦ 


rl5 

=  rl**5 

dO 

r25 

=  r2**5 

dO 

rl7 

=  rl**7 

dO 

r27 

=  r2**7 

do 

c 


go 

to 

(30001, 

30002, 

30003, 

30004), i 

30001 

go 

to 

(30110, 

30120, 

30130, 

30140), j 

30002 

go 

to 

(30210, 

30220, 

30230, 

30240), j 

30003 

go 

to 

(30310, 

30320, 

30330, 

30340) , j 

30004 

go 

to 

(30410, 

30420, 

30430, 

30440), j 

c  note 

matrix  is 

quite  sparse  now . 

30110 

go 

to 

(30111, 

30112, 

30113, 

30114),)?, 

30130 

go 

to 

(30131, 

30132, 

30133, 

30134), k 

30310 

go 

to 

(30311, 

30312, 

30313, 

30314), k 

30330 

go 

to 

(30331, 

30332, 

30333, 

30334), k 

30111  h  =  -9.d0*xmua*qa/rl5  -  9.d0*xmu*qb/r25 

1  +  15.d0*xmua*qa*qa*qa/rn  +  15.d0*xmu*qb*qb*qb/r27 

return 

30112  h  =  O.dO 
return 

30113  h  =  -3 .  d0*xmua*x  (3) /rl5  -  3  .d0*xn\u*x  (3) /r25 

1  +  15.d0*xmua*qa*qa*x (3) /rll  +  15.d0»xmu*qb*qb*x {3) /r27 

return 

30114  h  =  O.dO 
return 


c 

30120  h  =  O.dO 
return 

c 

30131  go  to  30113 

30132  h  =  O.dO 
return 

30133  h  =  -3  .d0*xmua*‘qa/rl5  -  3.d0*xinu*qb/r25 

1  +  15.d0*xmua*qa*x (3) *x (3) /rl7  +  15.d0»xmu*qb*x (3) *x (3) /r27 

return 
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30134  h  =  O.dO 
return 
c 

30140  h  =  O.dO 
return 

30210  h  =  O.dO 
return 

30220  h  =  O.dO 
return 

30230  h  =  O.dO 
return 

30240  h  =  O.dO 
return 
c 

30311  go  to  30113 

30312  h  =  O.dO 
return 

30313  go  to  30133 

30314  h  =  O.dO 
return 

c 

30320  h  =  O.dO 
return 
c 

30331  go  to  30133 

30332  h  =  O.dO 
return 

30333  h  =  -9.d0*xmua*x (3) /rl5  -  9.d0*xmu*x (3) /r25 

1  +  IS.dO* (xmua/rl7  +  xmu/r27) *x (3) *x (3) *x (3) 

return 
c 

30334  h  =  O.dO 
return 

30340  h  =  O.dO 
return 

30410  h  =  O.dO 
return 

30420  h  =  O.dO 
return 

30430  h  =  O.dO 
return 

30440  h  =  O.dO 
return 
c 

end 
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c 

subroutine  leqt2f (a, m, n, nn, b, Idgt, x, ler) 
c 

c  gaussian  elimination  with  maximal  pivoting 

c  interface  simulates  IMSL  routine 

c  solution  of  a  system  of  linear  equations  for  m  right  sides 

c  a:  matrix  of  system 

c  m:  number  of  rhs 

c  n:  order  of  a,  rows  in  b 

c  ia:  row  dimension  of  a,b 

c  b:  right  hand  sides ....  solution  on  return 

c  idgt:  ignored  here.... in  imsl  0=no  acc  test  on  input 

c  ldgt=  Idigits  ok  on  output  in  imsl 

c  x:  in  imsl,  n**2  +  3*n 
c  ier:  129:  singular  matrix,  0=ok 

c 

dimension  a (nn, nn) , b (nn,m) , irr (50) , x (1) 
double  precision  a, b, x, anorm, amax, p, tol 
c 

c  find  max  norm  of  a 

anorm  =  O.dO 
do  5  i  =  1 ,  n 

do  5  j  =  1 ,  n 

if  (dabs  (a  (i,  j)  )  .gt.  anorm)  anorm  =  dabs(a(i,j)) 

5  continue 

c  set  tolerance  =  2**  (-  number  of  binary  digits  in  mantissa) 

tol  =  l.d-12 
ler  =  0 
id  =  1 

do  10  i  =  l,n 
Irr(i)  =  0 
10  continue 
20  ir  =  1 
is  =  1 
amax  =  O.dO 
c  find  max  pivot 
do  60  i  =  l,n 

lf(irr(i))  60,30,  60 
30  do  50  j  =  l,n 

p  =  dabs (a (i, j) ) 
if(p-amax)  50,50,40 
40  ir  =  i 

is  =  j 
amax  =  p 
50  continue 

60  continue 

c  singularity  test 

if (amax/anorm  .gt.  tol)  go  to  70 
ier  =  129 
go  to  120 

c  forward  elimination 

70  irr(ir)  =  is 
do  90  i  =  l,n 

if(i  .eq.  ir)  go  to  90 
p  =  a(i,is)/a(ir,is) 
do  80  j  =  l,n 

a(i,j)  =  a(i,j)  -  p*a(ir,j) 

80  continue 

a(i,ls)  =0.0 
do  85  j  =  l,m 

b(i,j)  =  b(i,j)  -  p*b{ir,  j) 

85  continue 

90  continue 
id  =  id  +  1 

if(ld  .le.  n)  go  to  20 
c  back  substitution 
do  115  j  =  l,m 

do  100  i  =  l,n 
ir  =  irr(i) 

x  (ir)  =  b(l,  j)  /a  (i,  ir) 

100  continue 

do  110  i  =  l,n 

b(i,j)  =  x(i) 

110  continue 

115  continue 
120  return 
end 
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c 

subroutine  cleqt2f (a, m, n, nn, h, Idgt, x, ier) 
c 

c  complex  version  of 

c  gaussian  elimination  with  maximal  pivoting 

c  interface  simulates  IMSL  routine 

c  solution  of  a  system  of  linear  equations  for  m  right  sides 

c  a:  matrix  of  system 

c  m:  number  of  rhs 

c  n:  order  of  a,  rows  in  b 

c  ia:  row  dimension  of  aib 

c  b:  right  hand  sides ....  solution  on  return 

c  idgt:  ignored  here.... in  imsl  0=no  acc  test  on  input 

c  idgt=  #dlgits  ok  on  output  in  imsl 

c  x:  in  imsl,  n**2  +  3*n 

c  ier:  129:  singular  matrix,  0=ok 

c 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 
complex* 16  a (nn, nn) , b (nn, m) , x (28) , check, d 
dimension  irr (50) , rcheck (2) 
equivalence  (check, rcheck) 
c 

c  find  max  norm  of  a 

anorm  =  O.dO 
do  5  i  =  1 ,  n 

do  5  j  =  1 ,  n 

check  =  a (1, j) 

temp  =  dsqrt (rcheck ( 1) *rcheck (1) 

1  +rcheck (2) *rcheck  (2) ) 

if (temp  .gt.  anorm)  anorm  =  temp 

5  continue 

c  set  tolerance  =  2**  (-  number  of  binary  digits  in  mantissa) 

tol  =  l.d-12 
ier  =  0 
id  =  1 

do  10  i  =  l,n 
irr(i)  =  0 
10  continue 
20  ir  =  1 
is  =  1 
amax  =  O.dO 
c  find  max  pivot 
do  60  1  =  l,n 

if(irr(i))  60,  30,  60 
30  do  50  j  =  l,n 

check  =  a  (i,  j) 

p  =  dsqrt (rcheck (1) *rcheck (1) 

1  +rcheck (2) *rcheck (2) ) 

if(p-amax)  50,50,40 
40  ir  =  i 

is  =  j 
amax  =  p 
50  continue 

60  continue 

c  singularity  test 

if (amax/anorm  .gt.  tol)  go  to  70 
ier  =  129 
go  to  120 

c  forward  elimination 

70  irr(ir)  =  is 
do  90  i  =  l,n 

if(l  .eq.  ir)  go  to  90 
d  =  a  (i,  is) /a  (ir,  is) 
do  80  j  =  l,n 

a(i,j)  =  a(i,j)  -  d*a(ir,j) 

80  continue 

a (i, is)  =  0.0 
do  85  j  =  l,m 

b(i,j)  =  b(i,j)  -  d*b(ir,j) 

85  continue 

90  continue 
id  =  id  +  1 

if(id  .le.  n)  go  to  20 
c  back  substitution 

do  115  j  =  l,m 

do  100  i  =  l,n 
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ir  =  i rr (i) 

x(ir)  =  b(i,j)/a(i,  ir) 
100  continue 

do  110  i  =  l,n 

b(i,  j)  =  x(i) 

110  continue 

115  Continue 
120  ret  irn 
end 
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c 

subroutine  check (phi, val, vec, period, error) 
c 

c  -  checks  the  transformation  from  phi  to  jordon  normal  form 
c 

implicit  double  precision  (a-h) 

Implicit  integer  (i-n) 
implicit  double  precision  (o-z) 
dimension  xww (2) , xwc (2) , phi (4, 4 ) 
complex* 16  ww,  wc,  e jt (4, 4) , vec (4, 4 ), val  (4) 
complex*  16  fe(4,4),pf(4,4),pi(4,4),piv(4,4) 
equivalence  (ww,xww) 
equivalence  (wc,xwc) 
c 

c  *****************  calculate  e^ (Jt)  ********************** 

c 

do  100  i  =  1,4 

do  120  j  =  1,4 

ejt(j,l)  =  (O.dO,  O.dO) 

120  continue 

ejt  (i, i)  =  val (i) 

100  continue 

if(ejt(3,3)  .eq.  (l.dO,  O.dO))  then 
xww(l)  =  period 
xww (2)  =  O.dO 
ejt(3,4)  =  ww 

else 

endif 

c 

c  **************  FIND  phi*f  -  f*e''f.,-.)  ********************** 

c 

do  200  i  =  1, 4 

do  200  ]  =  1,  4 

fe(i,j)  =  (O.dO,  O.dO) 
pf  (i,  j)  =  (O.dO,  O.dO) 
do  200  k  =  1,  4 

fe(i,j)  =  fe(i,j)  +  vec (i, k)  *e jt  (k,  j) 
pf(i,j)  =  pf(i,j)  +  phi  (i,k)  *vec(k,  j) 

200  continue 

error  =  (O.dO,  O.dO) 
ww  =  error 
do  220  i  =  1,4 

do  220  j  =  1, 4 

wo  =  pf  ( j,  i)  -  fe  ( j,  i) 

if  (dabs  (xwc  (1) )  .gt.  dabs  (xww(l) ) )  xww(l)  =  xwc(l) 
if (aabs (xwc (2) )  .gt.  dabs (xww (2) ) )  xww(2)  =  xwc(2) 

220  continue 
c 

C  ******************  find  (phi-e-' (Jt)  *I)  *vec  **************** 
c 

do  300  i  =  1,4 

do  320  j  =  1,4 

do  340  k  =  1,4 

pi (k, j)  =  phi (k, j) 
piv(j,k)  =  (O.dO,  O.dO) 

340  continue 

pi(j;j)  =  phi(j,j)  -  val(i) 

320  continue 

do  360  j  =1,4 

do  360  k  =1,4 

piv(j,i)  =  piv(j,i)  +  pi  ( j,k)  *vec(k,  i  ) 

360  continue 

300  continue 

if(ejt(3,4)  .eq.  (O.dO,  O.dO))  goto  400 

piv(2,4)  =  piv(2,4)  -  vec  (2,  3)  *period 
piv(3,4)  =  piv(3,4)  -  vec  (3,  3)  *period 
do  400  i  =  1,4 

do  400  j  --1,4 

wc  =  piv(j,i) 

if  (dabs  (xwc  (1)  )  .gt.  dabs  (xww(l) ) )  xww(l)  =  xwc(l) 
if (dabs (xwc (2) )  .gt.  dabs (xww (2) ) )  xww(2)  =  xwc(2) 

400  continue 
error  =  ww 

2  format (lx, il, 2 (2x,  d24 . 17) ) 

3  format (2x, 2 (2x, d24 . 17) ) 
return 

end 
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c 

c 

c 

c 


c 

c 


c 


c 


c 

c 

c 


subroutine  symplec  (vec,  tzvec,  e  'ror) 

-  checks  if  eigenvector  matrix  is  syraplectic 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 


dimension  z (4, 4) , xww (2) , xwc (2) 

complex*  16  vec  (4,  4) ,  tvec  (4,4),  zvec  (4,  4) ,  tzvec  (4,  4) 
complex*16  error, ww,wc 

equivalence (ww, xww) 
equivalence (wc, xwc) 


data  z/ 

O.dO, 

-1  .dO, 

O.dO, 

O.dO, 

1 

1  .dO, 

O.dO, 

O.dO, 

O.dO, 

2 

O.dO, 

O.dO, 

O.dO, 

-l.dO, 

3 

O.dO, 

O.dO, 

l.dO, 

O.dO/ 

********************  transpose  vec. 


STORE  AS  tvec  ******************** 


do  440  i  =  1, 4 

do  440  j  =  1,4 

tvec(i,  j)  =  vec(j,  i) 

440  continue 


c 

c  ***************  CALCULATE  Z  VEC  ********************* 
c 

do  480  i  =  1,4 

do  480  j  =  1, 4 

zvec(i,j)  =  (O.dO,  O.dO) 
do  480  k  =  1,4 

zvec(i,j)  =  zvec(i,3)  +  z  (i, k)  *vec (k,  j) 

480  continue 
c 

c  T 

c  **************  CALCULATE  VEC  *  Z  VEC  =  Z  ************* 
c 

do  500  i  =  1,4 

do  500  j  =  1,4 

tzvec{i,j)  =  (O.dO,  O.dO) 
do  500  k  =  1,4 

tzvec(i,j)  =  tzvec(i,j)  +  tvec  (i,  k)  *zvec  (k,  j) 

500  continue 
c 

c  calculate  max  error 
c 

error  =  (O.dO,  O.dO) 
ww  =  error 
do  600  i  =  1,4 

do  600  j  =  1,4 

wc  =  tzvec(j,i)  -  z(j,i) 

if  (dabs  (xwc  (1)  )  .gt.  dabs  (xww(l) )  )  xww(l)  =  xwc(l) 
if (dabs (xwc (2) )  .gt.  dabs (xww (2) ) )  xww (2)  =  xwc (2) 

600  continue 
error  =  ww 


2  format  (lx,  i  1,  2  (2x,  d24 . 17)  ) 

3  format  (2x,  2  (2x,  d24 . 17)  ) 


return 

end 
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c 

c 

c 

c 


c 

c 


c 


c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 


subroutine  real (vec, x j, vecnew, x jnew, itype) 

-  converts  imaginary  E  and  J  to  real  E  and  J 

implicit  double  precision  (a-h) 
implicit  integer  (1-n) 
implicit  double  precision  (o-z) 

dimension  xwc (2) , xww (2) , itype (2) 

complex*  16  vec(4,4),xj(4,4),  vecnew  <4,  4) ,  x  jnew  (4,4) 
complex*  16  ww,wc,t(4,4),ti(4,4),  temp  (4,  4) 
complex*16  vec2 (4, 4) , x j2 (4, 4) 

equivalence  (wc,xwc) 
equivalence  (ww, xww) 

*************  remove  imaginary  POINCARE  EXPONENTS  ******** 

****************  CALCULATE  T  AND  T  INVERSE  **************** 

do  20  i  =  1,4 

do  25  j  =  1,4 

t(i,j)  =  (O.dO,  O.dO) 
ti(i,j)  =  (O.dO,  O.dO) 

25  continue 

t(i,i)  =  (l.dO,  O.dO) 
ti(i,i)  =  (l.dO,  O.dO) 

20  continue 

do  40  i  =  1,3,2 

iii  =  int ( (1+1) /2) 
if  (itype(iii)  .eq.  2)  then 
t(l,i)  =  (.5d0,  -.5d0) 
t(i,i  +  l)  =  (.5d0,  -.5d0) 
t(i  +  l,i)  =  (-.5d0,  -.5d0) 
t(i+l,i  +  l)  =  (.5d0,  .5d0) 
ti(i,i)  =  (.5d0,  .5d0) 
ti(i,i  +  l)  =  (-.5d0,  .5d0) 
ti(i  +  l,i)  =  (.5d0,  .5d0) 
ti(i  +  l,i  +  l)  =  (.5d0,  -.5d0) 

else 

endif 

40  continue 

**************  CALCULATE  NEW  E  AND  J  MATRICES  ************ 

do  60  i  =  1,4 

do  60  j  =  1,4 

temp(i,j)  =  (O.dO,  O.dO) 
do  60  )^  =  1,4 

temp(i,j)  =  temp(i,j)  +  t  (i,  k)  *x  j  (X,  j) 

60  continue 

do  80  i  =  1,4 

do  80  j  =  1,4 

vec2(i,j)  =  (O.dO,  O.dO) 
xj2(l,j)  =  (O.dO,  O.dO) 
do  80  k  =  1,4 

vec2(i,j)  =  vec2(i,j)  +  vec  (i,  k)  *ti  (k,  j) 
xj2(i,j)  =  xj2(i,j)  +  temp(i,  k)  *11  (k,  j) 

80  continue 

****************  remove  imaginary  eigenvectors  **************** 

****************  CALCULATE  T  AND  T  INVERSE  **************** 

do  100  i  =  1,4 

do  105  j  =  1,4 

t(i,j)  =  (O.dO,  O.dO) 
ti(i,j)  =  (O.dO,  O.dO) 

105  continue 

t(i,i)  =  (l.dO,  O.dO) 
ti(i,i)  =  (l.dO,  O.dO) 

100  continue 

do  120  i  =  1,3,2 
itest  =  0 
do  140  j  =  1,4 

wc  =  vec2(j,i) 


140 


if  (dabs (xwc (2) )  -gt.  l.d-10)  itest  =  1 

140  continue 

if  (itest  .eq.  1)  then 

t(i,i)  =  (O.dO,  O.dO) 
t(i,i  +  l)  =  (O.dO,  l.dO) 
t(i  +  l,i)  =  (O.dO,  l.dO) 
t(i  +  l,i  +  l)  =  (O.dO,  O.dO) 
ti(i,i)  =  (O.dO,  O.dO) 
ti(i,i+l)  =  (O.dO,  -l.dO) 
ti(i  +  l,i)  =  (O.dO,  -l.dO) 
ti(i  +  l,i  +  l)  =  (O.dO,  O.dO) 
iii  =  int ( (i+1) /2) 

if (itype (iii)  .eq.  0)  itype(iii)  =  -1 

else 

endif 

120  continue 
c 

c  ********  CALCULATE  NEW  EIGENVECTOR  MATRIX  AND  NEW  J  MATRIX  ****** 
C 

do  180  i  =  1,4 

do  180  j  =  1, 4 

temp(i,j)  =  (O.dO,  O.dO) 
do  180  )c  =  1,4 

temp(i,j)  =  temp(i,j)  +  t  (i,  k)  *x  j2  (k,  j) 

180  continue 

do  200  i  =  1,4 

do  200  j  =  1,4 

vecnew(i,j)  =  (O.dO,  O.dO) 
xjnew(i,j)  =  (O.dO,  O.dO) 
do  200  k  =  1,4 

vecnew(i,j)  =  vecnew(i,j)  +  vec2  (i,  k)  *ti  (k,  j) 
xjnew(i,j)  =  xjnew(i,j)  +  temp  (i,  k)  »ti  (k,  j ) 

200  continue 

return 

end 
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c 

c 

c 

c 

c 


c 

c 


c 

c  *** 

c 

c 

c  *** 
c 


200 

c 

c  *** 

c 

400 

500 

c 

c  *** 
c 


subroutine  fourier (f , ck, sk, n) 

subroutine  to  create  set  of  100  fourier  coefficients  for  a 
given  periodic  variable 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 

dimension  f (200) , ck (101) , sk (101) 

pi  =  dacos (-1 . dO) 
twopi  =  2.d0*pi 
alpha  =  twopi/dble (2*n) 
n2ml  =  2*n-l 

**********  ORDER  K  LOOP  ****************** 

do  500  k  =  0,n 

**********  COSINE  SUM  ****************** 

ck(k+l)  =  O.dO 
do  200  j  =  0,n2ml 

ck(k+l)  =  ck(k+l)  +  f ( j+1) *dcos (dble (k* j) *alpha) 
continue 

ck(k+l)  =  ck (k+1) /dble (n) 

*********  SINE  SUM  ************************ 

if(k  .eg.  0)  goto  500 
if(k  -eg.  n)  goto  500 
sk(k+l)  =  O.dO 
do  400  j  =  l,n2ml 

sk(k+l)  =  sk(k+l)  +  f < j+1) *dsin (dble (k*j) *alpha) 
continue 

sk(k+l)  =  sk (k+1) /dble (n) 
continue 

:*******  CORRECT  FIRST  AND  LAST  COSINE  COEFFICIENT  ****** 

ck(l)  =  .5d0*ck(l) 
ck(n+l)  =  .5d0*ck(n+l) 

return 

end 
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c 

c 

c 

c 

c 


c 

c 


c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 

c 


subroutine  rhs(k) 

canonical  EOM  and  EOV,  4th  order  system 
Phi  matrix  stored  by  cols,  end  to  end 

common  /data/  xmu,  xmua 

common  /ham/  t,  x  (20,  4) ,  f  (20,  4) ,  err  (20)  ,hh,  nn,mode 
external  h 

implicit  double  precision  (a-h) 
implicit  integer  (i-n) 
implicit  double  precision  (o-z) 

dimension  x  (20,  4) ,  f  (20,  4),err(20),xx(4),z(4,4)  ,grdh  (4,4),A(4,4) 


data  z/  O.dO, 

-l.dO, 

O.dO, 

O.dO, 

l.dO, 

O.dO, 

O.dO, 

O.dO, 

O.dO, 

O.dO, 

O.dO, 

-l.dO, 

O.dO, 

O.dO, 

l.dO, 

O.dO/ 

do  10  i  =  1,4 
10  xx(i)  =  x(i,k) 

******************  equations  of  motion  ******************** 

f(l,k)  =  h  (XX,  1,  2,  0,  0,  0,  0) 
f(2,k)  =  -h(xx,  1,  1,0, 0,0,0) 
f(3,k)  =  h  (xx,  1,  4,  0,  0,  0,  0) 
f(4,k)  =  -h(xx,l,3,0,0,0,0) 

if (mode  .eg.  0)  return 

***********  CALCULATE  ORDER  2  GRADIENT  MATRIX  ************* 

do  20  i  =  1,4 

do  20  j  =  1,4 

grdh(i,j)  =  h  (xx,  2,  i,  j,  0, 0,  0) 

20  continue 


c 

do  30  i  =  1,4 

do  30  ii  =  1,4 

A(i,ii)  =  O.dO 
do  30  j  =  1,4 

A(i,ii)  =  A(i,ii)  z  (i,  j)  *grdh  ( j,  ii ) 

30  continue 
c 

c  *********************  CALCULATE  A  PHI  ***************** 

c 

c  row  loop 

do  35  i  =  1,4 
c  col  loop 

do  35  ii  =  1,4 

ij  =  4*ii+l 
f(ij,k)  =  O.dO 
do  35  j  =  1,4 

f(ij,k)  =  f(ij,k)  t  A(i,  j)  *x  (4*ii  +  j,  k) 

35  continue 
return 
end 
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o  o 


c 

c 

c 

c 

c 


c 


subroutine  rhs(k) 

canonical  EOM  and  EOV,  4th  order  system 
Phi  matrix  stored  by  cols,  end  to  end 

common  /data/  xmu, xmua 

common  /ham/  t,  x  (20,  4) ,  f  (20,  4) , err  (20) , hh,  nn,mode 
common  / f dat /  x  j ( 4 , 4 ) 

external  h 


implicit  double  precision  (a-h) 

Implicit  integer  (i-n) 
implicit  double  precision  (o-z) 

dimension  x(20,4),f(20,4),err(20),xx(4),z(4,4) 
dimension  grdh  (4,4),a(4,4),  vec  (4,4) 
dimension  x j (4, 4) , fdot (4, 4) 
dimension  fl (4, 4) , f2  (4, 4) 


data  z/  O.dO, 

-l.dO, 

O.dO, 

O.dO, 

l.dO, 

O.dO, 

O.dO, 

O.dO, 

O.dO, 

O.dO, 

O.dO, 

-l.dO, 

O.dO, 

O.dO, 

l.dO, 

O.dO/ 

c 

do  10  i  =  1,4 

XX  (i)  =  x(i,)c) 

10  continue 

c  *******************  equations  of  motion  ************** 

f  (l,)t)  =  h(xx,l,2,0,0,0,0) 
f(2,)t)  =  -h  (XX,  1,1,  0,0,  0,0) 
f(3,)t)  =  h(xx,l,4,0,0,0,0) 
f(4,)?,)  =  -h(xx,l,3,0,0,0,0) 

if (mode  .eq.  0)  return 

*************  CALCULATE  ORDER  2  GRADIENT  MATRIX  »**♦***» 


do  20  i  =  1,4 

do  20  j  =  1, 4 

grdh(i,j)  =  h (xx, 2,  i,  j,  0,  0,  0) 

20  continue 
c 

do  30  i  =  1,4 

do  30  ii  =  1,4 

a(i,ii)  =  O.dO 
do  30  j  =  1,4 

a(i,ii)  =  a(i,ii)  +  z  (i,  j)  *grdh  ( j,  ii) 

30  continue 

c  *****************  CALCULATE  A*F  AND  F*J  ***************** 
c 

do  33  i  =  1,4 

do  33  j  =  1,4  ♦ 

vec(j,i)  =  x(i*4+j,)c) 

33  continue 
c  row  loop 

do  35  i  =  1,4 
c  col  loop 

do  35  j  =  1,4 

fl  (i, j)  =  O.dO 
f2(i,  j)  =  O.dO 
do  35  ii  =  1,4 

fl(i,j)  =  fl(i,j)  +  a  (i,  ii)  *vec(ii,  j) 
f2(i,j)  =  f2(i,j)  +  vec(i,  ii)  *xj  (ii,  j) 

35  continue 

do  36  i  =  1,4 

do  36  j  =  1,4 

fdot(i,j)  =  fl(i,  j)-f2(i,3) 

36  continue 

do  37  i  =  1,4 

do  37  j  =  1,4 

f(i*4  +  j,)t)  =  fdot(j,l) 

37  continue 


return 

end 
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c 

subroutine  rhs(k) 
c 

c  calculate  rhs  for  nearly-periodic  eom,  using 

c  expanded  hamiltonian 

c 

c  canonical  EOM  and  EOV,  4th  order  system 

c 

common  /datarhs/  wO, wl, w2, ck (20, 100) , sk (20, 100) , itype (2) 
common  /ham/  t ,  x  (20,  4) ,  f  (20,  4)  ,  err  (20) ,  hh,  nn,  mode 
c 

implicit  double  precision  (a-h) 
implicit  integer  (1-n) 
implicit  double  precision  (o-z) 
c 

dimension  x(20,4),f(20,4),err(20) 
dimension  coss (100),3lnn(100), itype (2) , w (2) 
dimension  ck(20, 100)  ,sk(20, 100) ,  c  (20) ,  b  (4) ,  bdot  (4) 

W (1)  =  Ml 
w(2)  =  w2 
c 

c  **************  generate  SIN(1  to  50  *w0)AND  COS(l  to  50  *  wO) 
c 

coss(l)  =  dcos(w0*t) 

sinn(l)  =  dsin(w0*t) 

coss(2)  =  2.d0*coss (1) *coss (1) -l.dO 

sinn(2)  =  2 .d0*sinn (1) *coss (1) 

do  100  i  =  3,100 

coss(i)  =  2 .d0*coss (i-1) *coss (1) -coss (i-2) 
sinn(i)  =  2 .d0*sinn (i-1 ) *coss (1) -sinn (i-2) 

100  continue 


c 

c  **********  reconstruct  periodic  function  coefficients  ****************** 

c 

do  200  1  =  1,20 

c  (i)  =  ok  (i,  1) 
do  200  j  =  1, 99 

c(i)  =  c(i)  +  ck (i, j+1) *coss ( j)  +  sk (i,  j  +  1) *sinn { j ) 

200  continue 


c 

c  **************  establish  modal  vectors  ************************ 


c 

do  240  i  =  1,4 

b(i)  =  x(i,k) 

240  continue 
c 

C  **************  CALCULATE  BDOT  AND  ESTABLISH  EOM  ****************** 


c 


bdot ( 1 ) 

& 

& 

& 

& 

bdot (2) 

& 

& 

& 

& 

bdot  (3) 

& 

& 

& 

& 

bdot (4) 

& 

& 

& 

& 


=  +l.d0*b(l)  *b(l)  *c(2) 

+  l.d0*b(l)  *b(3)  *c(6) 

+  3.d0*b{2)»b(2)*c(ll) 
+  2.d0*b(2)  *b{4)  *c(13) 
+  l.d0*b(3)*b(4)*c(15) 
=  -3.d0*b(l)  *b(l)  *c(l) 

-  2.d0*b(l)  *b(3)  *c(3) 

-  I.d0*b(2)*b(2)»c(5) 

-  I.d0*b(2)  *b(4)  *c(7) 

-  I.d0*b(3)  *b(4)  *c(9) 

=  l.d0*b(l)*b(l)*c(4) 

+  l.d0*b(l)  *b(3)  *c(9) 

+  l.d0*b(2)  *b(2)  *c(13) 
+  2.d0*b(2)  *b(4)  *c(16) 
+  2.d0*b(3)  *b(4)  *c(19) 
=  -l.d0*b(l)  *b(l)  *c(3) 

-  2.d0*b(l)  *b(3)  *c(8) 

-  I.d0*b(2)  *b(2)  *c(12) 

-  I.d0*b(2)  *b(4)  *c(15) 

-  2.d0*b(3)  *b(4)  *c(18) 


+  2.d0*b(l)*b(2)*c(5) 

+  l.d0*b(l)*b(4)»c(7) 

+  2.d0*b(2)  *b(3)  *c(12) 
+  l.d0*b(3)*b(3)*c(14) 
+  l.d0*b(4)  *b(4)  *c  (16) 

-  2.d0»b(l)*b(2)  *c(2) 

-  2.d0*b(l)*b(4)*c(4) 

-  I.d0*b(2)*b(3)*c(6) 

-  I.d0*b(3)*b(3)*c(8) 

-  I.d0*b(4)  *b(4)  *c  (10) 
+  l.d0*b(l)*b(2)*c(7) 

+  2.d0*b(l)  *b(4)  *c(10) 
+  l.d0»b(2)  *b(3)  *0(15) 
+  l.d0*b(3)  *b(3)  *c(18) 
+  3.d0*b(4)  »b(4)  *c  (20) 

-  I.d0*b(l)*b(2)*c(6) 

-  I.d0*b(l)*b(4)*c(9) 

-  2.d0»b(2)  *b(3)  *0(14) 

-  3.d0»b(3)*b(3)*c(17) 

-  I.d0*b(4)*b(4)*c(19) 


do  250  i  =  1,3,2 

iii  =  int  (  (i  +  1) /2) 
if  (Itype(iii)  .eq.  -1)  then 

bdot (i+1)  =  bdot (1+1)  +  b(i) 
elseif  (itype (iii)  .eq.  0)  then 
bdot(i)  =  bdot(i)  +  b(i+l) 
elseif  (itype(iii)  .eq.  1)  then 

bdot(i)  =  bdot(l)  +  w(lii)*b(i) 
bdot(i  +  l)  =bdot(i  +  l)  -  M(lil)  *b(l  +  l) 
elseif  (itype(lii)  .eq.  2)  then 
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bdot(i)  =  bdot(i)  +  w  (iii)  *b  (i+1) 
bdot(i+l)  =  bdot(l+l)  -  w(iii)*b(i) 

else 

endlf 

250  continue 

do  260  i  =  1,4 

f  (i,k)  =  bdot  (i) 

260  continue 

return 

end 
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